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CHAPTER  1 

INTRODUCTION 

Nickel-cadmium  storage  catteries  have  been  in  use  for 
more  than  fifty  years;  since  World  War  II,  r.e w  design  princi¬ 
ples  have  been  introduced  which  have  greatly  extended  their 
usefulness.  Many  theories  and  investigations  have  been  put 
forward  concerning  the  reaction  mechanisms  of  both  cadmium  and 
nickel  electrodes.  Because  the  activity  of  the  nickel  electrode 
has  been  more  fully  characterized  than  that  of  the  cadmium 
electrode,  we  will  focus  our  attention  on  the  behavior  of  the 
cadmium  electrode. 

There  are  a  number  of  methods  used  for  the  fabrication 

of  cadmium  electrodes.  These  methods  are  usually  classified 

into  three  categories!  (a)  pressed  powder  or  pasted  typev ^ , 

(4  5  6  7  S ) 

(b)  impregnated  nickel  sinter  typev  ,-3t  ’  '  '  and  (c)  electro- 

(9  10  11  12 ) 

chemical  impregnation  71  ’  ’  The  last  method,  however,  is 

simple  in  operation  and  can,  in  theory,  achieve  high  loading  of 
active  material  in  a  single  impregnation  cycle  at  the  lowest  cost. 

In  the  electrochemical  impregnation  process,  electro- 
chemically  active  material  is  introduced  into  the  porous  nickel 
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plaqe  by  electrochemical  precipitation.  The  plaque,  which  serves 
as  a  cathode,  is  submerged  into  an  electrolysis  cell  containing 
an  aqueous  solution  of  cadmium  nitrate  at  a  temperature  near 
the  boiling  point  of  the  solution.  The  initial  pH  value  of  the 
impregnation  solution  is  adjusted  to  around  3  by  titration  with 
nitric  acid.  (As  the  pH  value  of  the  bath  increases,  i.e., 
becomes  more  basic,  cadmium  hydroxide  begins  to  precipitate  in 
regions  other  than  within  the  pores  of  the  plaque,  and  active 
material  is  lost  in  the  bath,  which  is  not  what  we  want.) 

The  electrochemical  reaction  at  the  cathode  during  impre¬ 
gnation  could  either  be  the  liberation  of  hydrogen, 

2H£0  +  2e“  - ►  H2  +  20H"  ,  (1) 

or  the  reduction  of  nitrate  ion, 


(2) 


3oth  reactions  produce  hydroxy  ions  which  then  either  diffuse 
away  from  the  electrode  porous  surface  or  react  with  the  cadmium 
ions  present  in  the  pores  to  form  cadmium  hydroxide.  The  cadmium 


H20  +  NO^  +  2e“ - »  N02"  +20H 

•  •  •  •  • 

2H20  +  NO^"  +  2e”  -  HNC>2  +30H 

#  9  •  a  t 

6H20  +  ;i03"  +  8e“ - -  NH  +90H 
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hydroxide  then  deposits  within  the  pores  of  the  porous  plaque. 

The  precipitation  reaction  is: 

Cd+2  +  20H“  - *  Cd(  OH)  2 .  (3) 

There  is  the  possibility  of  forming  cadmium  according  to  the 
reduction  reaction: 

Cd(0H)2  +  2e“  - *-  Cd  +  20H".  (4) 

This  reaction  is  presumably  prevented  from  taking  place  by 
hydroxy  ion  generated  on  the  electrode  surface  by  reaction  (1) 
or  (2),  so  that  there  is  no  chance  for  cadmium  ions  to  diffuse 
through  the  solution  onto  the  surface  and  be  reduced  to  cadmium 
metal . 


Following  impregnation,  the  electrode  is  cathodically 
polarized  in  a  25  $  potassium  hydroxide  solution  with  the  impre¬ 
gnated  plate  connected  as  the  cathode  and  subjected  to  low  direct 

2 

anodic  current  (0.25  to  0.75  amp/in  ).  The  solution  is  maintained 
at  about  100°C  for  a  period  of  from  15  to  45  minutes  to  convert 
cadmium  hydroxide  to  cadmium  in  the  pores  of  the  plaque  by 
reaction  (4)  or  to  remove  any  nitrate  ion  or  other  undesirable 
reduction  products.  The  active  material  of  the  electrode  is  pre¬ 
sumed  to  be  both  cadmium  hydroxide  and  cadmium. 
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In  this  work,  the  electrochemical  impregnation  process 
for  the  fabrication  of  the  cadmium  electrode  was  examined  both 
experimentally  and  theoretically. 

Experimentally,  a  set  of  current-density- vs.-time  transients 
at  constant  applied  potential  were  obtained  for  a  nickel  micro¬ 
electrode.  The  current  density  was  due  to  the  occurrence  of  the 
electrode  reaction  (2)  and  precipitation  reaction  (3)  presented 
earlier.  For  the  purpose  of  identifying  the  electrochemical 
reactions,  another  set  of  current-density -vs.-time  transients  were 
obtained  using  an  electrolyte  which  contained  no  cadmium  ion. 

In  this  case,  the  precipitation  reaction  (3)  no  longer  took 
place.  The  reactions  at  the  electrode  were  the  only  reactions 
taking  place. 

The  transient  current  responsed  in  an  electrolyte  which 
contained  cadmium  ions  differed  significantly  from  that  in  an 
electrolyte  which  did  not  contain  the  cadmium  ions.  This  diff¬ 
erence  is  presumed  to  be  due  to  the  precipitation  reaction. 

The  current  obtained  in  the  presence  of  cadmium  ions  in  the 
electrolyte  at  any  time  was  higher  than  that  obtained  at  the 
same  applied  potential  in  the  absence  of  cadmium  ions.  This 
result  is  explainable  by  the  fact  that  the  cadmium  ions  in  the 
solution  consumed  the  hydroxy  ions  which  were  produced  by  the 
heterogeneous  reaction  (2)  and,  as  a  result,  promoted  the  hetero¬ 
geneous  electrode  reactions. 
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The  experimental  current-density-  vs.-time  curves  obtained 
in  the  potential-step  experiment  were  used  to  identify  both 
the  heterogeneous  electrochemical  and  the  homogeneous  preci¬ 
pitation  reaction  rate  expressions.  The  heterogeneous  rate  was 
identified  first  from  the  experimental  results  in  which  there 
was  no  cadmium  ion  in  the  electrolyte.  The  precipitation  rate 
was  then  identified  from  the  experimental  results  in  which  the 
cadmium  ions  were  present  in  the  electrolyte. 

A  set  of  differential  equations  which  described  the  trans¬ 
port  process  in  the  cadmium-ion-free  solution  was  solved  in  terms 
of  two  unknown  constants.  The  unknown  constants,  related  to  the 
heterogeneous  electrode  reactions,  were  identified  by  matching 
the  theoretical  results  with  the  results  of  the  experiment.  The 
precipitation  reaction  rate  constant  was  then  obtained  by  a 
similar  procedure. 
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CHAPTER  2 


LITERATURE  REVIEW 

Impregnation  of  nickel  sinters  with  cadmium  hydroxide, 
followed  by  electrochemical  reduction  to  cadmium,  is  the  most 
economical  method  by  which  to  fabricate  cadmium  electrodes. 

Almost  all  of  the  research  works  concerning  this  subject  may 
be  divided  into  three  types:  (a)  manufacturing  techniques  de¬ 
velopment,  (b)  impregnation  process  study  arid  (c)  charge-discharge 
study . 


During  recent  years,  manufacturing  techniques,  especially 
those  involving  electrochemical  impregnation,  have  employed 
widely  varying  operating  conditions^ .  There  are 
several  parameters  that  affect  the  electrochemical  impregnation 
process  such  as: 

1)  Temperature:  Beauchamp  emphasized  the  fact  that  high 

impregnation  bath  temperature  keeps  the  size  of  cadmium 
hydroxide  crystals  small,  and  thereby  obtaining  a  more  active 
electrode . 

2)  pH  value:  In  Beauchamp’s  process  ,  the  pH  value  was 

closely  controlled  by  buffering  the  cadmium  nitrate  solution 

/  A  '  \ 

with  sodium  nitrite.  Pickett  °'  in  a  similar  process, 
maintained  the  pH  value  of  the  cadmium  nitrate  bath  between 
3-5  by  titrating  the  solution  with  proper  amount;  of  nitric 
acid , 
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3)  Current  density:  The  current  density,  based  on  apparent 

( 14) 

electrode  surface,  is  another  important  variable.  Beauchamp 
used  D.  C.  current  of  4-8  amps/in^  while  Bulan  ^-5)  used 

2 

alternating  current  of  1.2-1. 6  amps/in  and  claims  that  the 
electrode  impregnated  by  alterning  current  technique  retained 
more  of  its  capacity  after  cycling  as  compared  to  electrodes 
fabricated  by  other  methods.  Pickett presents  another 
technique  similar  to  the  above  and  claims  that  his  technique 
yields  higher  loading  cadmium  electrodes  and  achieves  about 
85%  utilization  of  the  active  material. 

The  electrochemistry  of  the  Cd/Cd(0H)2  electrode  in  con¬ 
centrated  alkaline  solution,  typically  KOH  as  in  the  Ni-Ca 
battery,  has  been  studied  extensively  both  theoretically  and 
experimentally  ^ ^ .  Bennion,  et  al.^^  in  their 

theoretical  development,  present  two  models,  the  solution- 
diffusion  model  and  film  model,  in  considering  the  failure 
mechanism  in  the  performance  of  secondary  batteries  using  metal/ 
metal  salt  couples.  They  conclude  that  one  mode  of  failure  of 
the  metal/metal  oxide  porous  electrode  is  the  blocking  of  the 

pore  or  the  complete  covering  of  surface  by  product  crystallities . 

(18) 

Falk'  obtained  X-ray  diffraction  patterns  from  electrodes 
submerged  in  electrolyte  during  charge  and  discharge  by  means 
of  a  special  test  cell.  He  claims  that,  during  charge,  Cd(CH)2 

is  successfully  transformed  into  Ca  metal  in  the  cadmium  electrode, 
but  even  after  a  strong  overcharge,  this  transformation  is  not 
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complete,  as  indicated  by  the  X-ray  pattern  at  this  stage. 

Furthermore,  Cd(OH),  is  the  end  product  during  discharge  and 

there  is  no  evidence  of  the  presence  of  a  CdO  film  which  is 

believed  to  be  responsible  for  the  passivation  of  the  electrode. 

( 19 ) 

Gross  and  Glocking  summarized  the  results  of  recent  research 
by  characterizing  the  negative  cadmium  electrode  in  a  nickel- 
cadmium  cell  and  evaluating  some  of  the  published  results  about 
which  researchers  still  disagree  among  themselves,  such  as  the 
mechanism  for  passivation,  which  is  an  important  phenomena  that 
affects  the  capacity  of  the  nickel-cadmium  cell. 

In  comparison  to  the  extensive  studies  on  the  electro¬ 
chemistry  of  the  Cd/Cd(OH)2  electrode  in  KOH  solution,  there 
is  little  information  available  concerning  the  electrochemical 
reactions  taking  place  during  the  electrochemical  impregnation 
process  which  uses  cadmium  nitrate  solution  as  electrolyte. 

The  only  research  works  known  to  the  author  are  the  electro- 

(21  22 ) 

cnemicai  impregnation  studies  using  cycling  voltammetry  ’ 

2  2 
They  used  a  0.013  cm  nickel  microelectrode  and  a  0.064  cm 

cadmium  microelectrode.  Some  important  conclusions  of  their 

studies  are  summarized  below: 

(1)  There  was  a  large  reduction  wave  at  -O.63V  vs  SGE( Saturated 
Calomel  Electrode)  which  appeared  only  during  the  initial 
negative  scan.  This  is  believed  to  be  the  reduction  of  a 

solution  species  to  liberate  0H~  ion  in  the  presence  of 

+2 

Cd  to  precipitate  cadmium  hydroxide  on  the  electrode 


surface . 
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(2)  A  second  reduction  process  was  observed  at  -0.62V  vs  SGE 
during  the  reverse  positive  scan.  This  is  due  to  the  reductive 
formation  of  metallic  cadmium  from  cadmium  oxide  or  hydroxide 
that  was  deposited  on  the  electrode  surface  during  the  pre¬ 
vious  forward  scan. 

(3)  Some  processes  (or  one  of  the  processes)  during  the  first 
cycle  effectively  passivated  the  electrode  to  further  electro¬ 
activity.  The  formation  of  the  gray  material  at  -O.63V  vs 
SCE,  discussed  in  (1),  was  the  most  probable  reason  for  the 
passivation  of  the  electrode. 
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CHAPTER  3 

THEORETICAL  APPROACH 

3 . 1  General 

In  voltammetry  at  constant  voltage,  the  current  through 
the  electrolytic  cell  is  measured  as  a  function  of  time  while 
the  potential  of  the  polarized  electrode  is  held  at  a  constant 
value.  The  current  is  affected  by  both  the  transfer  process  and 
the  electrochemical  and/or  chemical  reactions. 

A  plane  electrode  is  immersed  in  a  solution  which  contains 
a  substance  0.  At  the  imposed  potential,  0  is  reduced  to  another 
substance  R  at  the  electrode  surface,  according  to  the  following 
electrochemical  reaction: 

0  +  ne"  - ►  R  (5) 

The  current-density- vs.-time  curve  obtained  at  a  constant 
voltage  depends  on  both  the  kinetics  of  the  electrochemical 
reaction  and  the  rate  of  mass  transfer. 

Under  our  experimental  conditions,  diffusion  is  the  only 
mode  of  mass  transfer.  The  other  two  modes  of  mass  transfer, 
namely,  migration  and  convection,  are  insignificant.  Migration 

I 
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is  eliminated  by  the  presence  of  a  large  excess  of  supporting 
electrolyte  which  reduces  the  potential  field.  Convection  is 
avoided  due  to  the  conditions  that  (1)  the  solution  is  not 
stirred  and  (2)  the  duration  of  electrolysis  is  short,  so  that 
the  density  change  has  not  yet  become  a  significant  factor  and 
natural  convection  has  not  taken  place. 


The  dimension  of  the  cell,  as  compared  to  the  microelectrode 
used,  is  so  large  that  the  solution  may  be  regarded  as  extending 
to  infinity;  i.e.,  the  diffusion  process  is  semi-infinite. 
According  to  Fick's  first  law,  the  flux  of  the  substance  0  at 
a  distance  x  from  the  electrode  and  in  the  direction  perpen¬ 
dicular  to  the  electrode  is: 

^C0(x,t) 

N0  (x.t)  =  -  D0-“  (6) 

ox  • 

Conservation  of  the  species  0  leads  to  the  following  differential 
equation : 

?>Cn(x,t)  =  D  ^2c0(x,t) 

t  0  7)  x2  (?) 


where  Cq  is  the  concentration  of  substance  0,  which  is  a  function 
of  x  and  t.  Dq  is  the  diffusion  coefficient  of  substance  0, 
and  is  assumed  to  be  a  constant. 
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When  the  substance  0  is  reduced  at  a  very  high  rate,  by 
imposing  a  sufficiently  large  potential  step,  such  that  the 
concentration  of  this  substance  at  the  electrode  surface  is 
reduced  to  zero  soon  after  the  start  of  electrolysis ,  the 
boundary  condition  at  the  electrode/electrolyte  interface  for 
equation  (7)  is  Cg(0,t)=0  for  t>0.  Initially  before  the  electro¬ 
lysis,  the  solution  is  homogeneous,  and  the  concentration  Cq 
at  t=0  is  uniform  with  the  bulk  value.  The  other  boundary 
condition  is  that  C0  approaches  its  bulk  value  when  x  is 
sufficiently  large. 

The  electrolysis  current  is  equal  to  the  product  of  the 
charge  involved  in  the  reduction  of  one  mole  of  substance  0  by 
the  flux  of  this  substance  at  the  electrode  surface.  Thus  the 
current  density  (current  per  unit  electrode  surface  area)  is: 

1=  nFN  0 ( 0 , t ) ,  (8) 

where  n  is  the  number  of  electrons  involved  in  the  reduction 
reaction  .of  substance  0,  F  is  the  Faraday  constant,  and  NQ(0,t) 
is  the  flux  of  substance  0  for  x=0. 

The  solution  of  the  above  problem  is^^  : 

1=  nFD  *  cj  - 1-J— - -  (9) 

U  U  t  2 

f 

where  Cq  is  the  bulk  value  of  substance  0. 


13 


The  current  density  represented  by  equation  (9)  is  called 
the  limiting  current  density  of  electrochemical  reaction  (5)« 

This  is  the  maximum  current  density  which  can  be  achieved  by 
this  reaction.  It  was  found  that  if  nitrate  ion  is  the  only 
reducible  substance,  the  current  density  predicted  by  equation 
(9)  is  considerably  lower  than  the  experimental  value.  This 
contradition  is  resolved  if  one  realizes  that  under  the  experi¬ 
mental  conditions  water  is  also  reduced. 

At  any  other  less-negative  potential  the  current  will 
depend  on  the  electrode  kinetics.  Now  the  condition  thax  the 
reducible  substance  concentration  becomes  zero  after  the  appli¬ 
cation  of  the  potential  step  has  to  be  replaced  by  a  kinetics 
equation,  which  is  a  function  of  potential  and  the  concentrations 
of  surface  reactant  and  product.  No  general  solution  is  available 

except  for  the  case  when  the  reaction  rate  is  linearly  dependent 

,  .  (22) 
on  tne  concentrations  of  the  surface  reactant  and  product  J  . 

3.2  Description  of  the  Model 

The  real  porous  electrode  is  formed  by  sintering  nickel 
powder  on  a  nickel  grid  to  form  a  plaque.  The  plaque  serves  as 
the  cathode,  which  is  positioned  in  a  cadmium  nitrate  solution 
with  a  pH  value  between  3  to  5*  The  passage  of  the  current 
produces  hydroxy  ions  which  cause  the  cadmium  to  be  deposited 
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into  the  porous  pores.  The  transport  phenomena  associated 
with  electrochemical  reactions  in  a  porous  medium  is  very 
complicated.  The  approach  here  is  to  study  the  intrinsic 
reactions  which  take  place  inside  the  pores  in  a  simple 
flat  electrode.  The  plate  electrode  is  submerged  in  a 
semi-infinite  pool  of  electrolyte  so  that  the  mass  transfer 
problem  can  be  treated  as  a  one-dimensional  problem  in 
the  x-direction  (perpendicular  to  the  electrode  surface). 

In  the  solution  side,  there  is  a  thin  double-charge  layer 
near  the  electrode.  The  thickness  of  the  double  layer  is 

O 

on  the  order  of  10  to  100  A.  It  is  too  thin  to  be  probed 
adequately,  and  the  theory  of  the  diffuse  layer  is  a 
microscopic,  rather  than  a  macroscopic,  phenomena.  The 
mass  transfer  region  to  be  considered  here  is  thus  located 
outside  of  this  double  layer. 

The  electrodeposition  process  operates  on  the  principle 
that  the  hydroxy  ion,  which  is  generated  by  the  reduction 
process,  coprecipitates  with  the  cadmium  ion  in  the  solution 
to  form  cadmium  hydroxide  crystalloid.  The  approach  taken 
here  to  unravel  the  complicated  sequence  of  reactions  is 
to  form  a  model  which  describes  the  transport  processes 
and  reactions.  The  model  will  first  be  used  to  analyze 
the  experimental  results  obtained  under  the  condition  of 
no  pricipitation  reaction,  by  using  an  electrolytic  solution 
which  is  free  of  cadmium  ions.  The  heterogeneous  reduction 
reaction  will  be  determined.  The  transport  and  reacxion 


15 


problems  will  then  be  examined  in  the  case  when  the  solution 
contains  cadmium  ions.  The  rate  of  the  precipixation 
reaction  will  then  be  determined  by  matching  the  model's 
prediction  and  the  experimental  results. 


3 . 3  Formulation  of  Transport  Process  Equations 


Three  ionic  species  to  be  considered  in  the  electro- 

-  -  +2  . 

chemical  system  are  N0^»  OH  and  Cd  ions.  In  the  absence 
of  migration  and  convection,  the  flux  of  each  species  is 
expressed  as^^  s 


"i  -  -  Di 


,  1  =  1, <2, 3  , 


where 


is  the  flux  of  species  i  (moles/cm  -sec ) , 
is  the  diffusion  coefficient  of  species  i  (cm  /sec), 
Ch  is  the  concentration  of  species  i  (moles/cnr ) , 

which  is  a  function  of  both  x  and  t,  and  i=l,2,3 

-  -  +2 

corresponding  to  N 0 ^ ,  CK  and  Cd  ions. 

(24) 

The  equation  of  continuity  of  each  species  is 


?Ci  = 

"St 


d  Ni 


i=l,2,3  , 


where  S ^  is  the  rate  of  production  of  species  i 
( raoles/cm  -sec )  by  precipitation  reaction,  which  occurs 
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only  in  the  presence  of  cadmium  ions. 


Substitution  of  equation  (10)  into  (11)  gives 


?>t 


2 

3  c. 


=  Di' 


2):. 


2 


i=l,2,3 


(12) 


It  is  assumed  that  the  precipitation  reaction  (3) 
can  be  treated  as  a  homogeneous  reaction  and  is  linearly 

proportional  to  the  degree  of  supersaturation.  Therefore, 

the  consumption  rate  of  cadmium  ion  due  to  this  reaction 

is  ; 


R 


3 


(13) 


where  k  is  the  homogeneous  reaction  rate  constant  (sec-^), 
K sp  is  the  solubility  constant  of  Cd(0H)2 

and  has  the  value  of  2xl0-^  (moles^/cm^)  ^5) . 

3y  the  stoichiometric  relation  of  equation  (3)>  the  con¬ 
sumption  rate  of  hydroxy  ion  due  to  the  same  reaction  is  : 


=  2  R3  - 


-2k 


(c2)‘ 


(14) 


Substitution  of  equations  (13)  ana  (14)  into  equation  (12) 
leads  to  the  following  two  differential  equations  s 


7)t 


->c. 


~d'^2  =  2  k 


Ksp 


(c2)‘ 


(  15) 


2  C 

TC, 

D- - k  C- 

3}x2  I  3 


/  p  \ 
v  2; 


For  NO.  ion, 


R:  =  o. 


Equation  (12)  for  NO^  ion  may  be  reduced  to 


7>^1  2Tci 

- —  =  n  - — 

■at  i  \  2 

O  X 


Only  two  of  the  three  differential  equations  need 

to  be  solved  for  the  concentration  profiles,  due  to  the 

(24) 

electroneutrality  condition  s 


2C3-C1-C2=0 


Substitution  of  equation  (18)  into  equation  (15) 


gives  s 


~2>  C2 
it 


??C2 


Equations  (1?)  and  (19)  are  solved  simultaneously 

to  obtain  the  concentration  profiles  of  ar.d  C~,  i.e., 

-  .  +2 

i.'O^  and  OH  ions.  The  profile  can  then  be  obtained  by 
equation  (18)  from  the  known  profile^  of  and  CU . 


A  set  of  differential  equations,  which  are  simpler 
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than  equations  (17)  and  (19).  ars  obtained  when  the  elec 
trolyte  contains  no  cadmium  ion.  In  the  absence  of  cad¬ 
mium  ion,  the  precipitation  reaction  (3)  can  not  occur; 
thus,  FL  terms  are  all  equal  to  zero.  Equation  (12)  is 
then  reduced  to i 


Ih. 

at 


i=l »2  , 


(20) 


in  which  i=l,2  corresponds 


to  NG^ 


and  OH 


ions . 


Equation  (20)  is  written  for  and  C^s 


7><>i 

7>t 


= 


(21) 


(22) 


The  other  condition,  the  concentrations  of  various 
ions  at  the  electrode/electrolyte  interface  after  the 
start  of  the  current  flow,  has  to  be  described  by  the 
electrode  kinetics  of  the  reduction  reaction  of  N0“  ions 


3.^  Electrode  Kinetics 


The  possible  reduction  reactions  are  numerous,  as 
shown  in  equation  (2).  Since  in  e*..ch  reaction,  OH"  ion 
is  generated,  we  will  represent  the  production  of  OH" 
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ions  as  the  following  electrochemical  reaction : 

2H20  +  NO-j  t  2e“ - ►  HNOg  +  30H".  (23) 

Because  there  are  two  electrons  involved  in  the  reaction, 
it  is  likely  not  an  elementary  step.  Let's  assume  the 
reaction  comprises  two  elementary  steps,  each  step  involving 
the  transfer  of  one  electron,  and  a  dissociation  reaction: 


(slow) ,  (24a) 

(fast),  (24b) 

(fast),  (25) 


where  and  Kc2  are  cathodic  reaction  constants,  and  K  ^ 
and  Ka2  are  anodic  reaction  constants,  respectively.  is 
the  dissociation  reaction  constant.  (OH) 2  ion  is  an  unstable 
intermediate  ion  which  immediately  enters  the  next  reaction  (24b/ 
The  rates  of  these  elementary  steps  should  always  be  propor¬ 
tional  to  each  other.  For  the  elementary  steps  listed  above, 
reaction  (24b)  should  occur  or.ce  every  time  reaction  (24a)  occurs 


Let  I.  and  I,  denote  the  cathodic  current  densities 
of  reactions  (24a)  and  (24b),  and  r,  and  r2  denote  the 
reaction  rates  of  reactions  (24a)  and  (24'c),  respectively. 
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Then,  tne  relationship  between  r.  ,  I-  ar.c.  surface  ion 


concentrations  of  H  C  ~ ,  (CH)^  ,  which  app 

( 2  ^ ) 

step  (24a),  may  be  written  as  ~  : 


ear  m  e-emer.tar: 


-  _Ll 

‘l“  nF 


(n=l ; 


=  Kal  (CN°c;](C?0H)-)  fJ 

W  W 


[l-P,  )? 


..  r„o  1  f  o  *!  r  l 

ci  "ho:  c:-:9o  e'‘?  "  ~ 

L  JJ  4  ^  *'  L  *" 


r  a?  „ 


Similarly,  for  reaction  (24b),  one  has: 


(n=l ) 


^K-r«pp^v] 

Ko2  [C(0K)-]  ex?[-^Lvj' 


where  p,  and  ^  are  symmetry  factors  representing  the  fractions 
of  applied  potential  V,  which  promotes  the  cathodic  reaction. 
C.°.~-  ,  and  C~u-  are  the  surface  ion  concentrations  of 

HOy  ( OH ) ^  and  OH”  on  the  electrode  surface. 


Since  reactions  (24a)  and  (2^c)  occur  at  the  same  rate, 
we  have : 


1=  1^  t  I ^  -  2i2  =  ’ 


2 
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where  I  is  overall  current  density.  It  it  further  assumed 
that  reaction  (24b)  is  fast,  and  is  essentially  in  equilibrium. 
From  equation  (27), 


Xa2  -y 

Xc2  exP 


and  = 
r 


=r  =  overall  reaction  rate. 


Substituting  equation  (28)  into  equation  (25)  and 
combining  equation  (29)  gives: 


-  KalKa2  (  ro 


[C°O-](C0H-]2  exp[' 


(2-A;f 


^  (c^J[cv]exp  [■  . 


From  reaction  (25)  one  may  write 


Xd  = 


(GHN02][C0H~] 

[GN0“][CH2gJ 


Substitution  cf  the  above  expression  into  equation  (jC)  leads 


22 
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ions  and  electrolyte  containing  no  cadmium  ions,  tv; o 
and  four  boundary  conditions  are  needed  to  solve  the 
equations  (17)  and  (19)  or  (21)  and  (22). 


initial 
differ en 


At  t=0,  the  concentration  of  each  ion  is  equal  to  its 
bulk  value  at  each  point  along  the  x-direction.  The  initial 
conditions  are  as  follows: 

1.  at  t=0,  C^=  for  all  x, 

2.  at  t=0,  C2=  C2  for  all  x, 

where  the  superscript  "c"  denotes  the  bulk  condition. 


The  first  boundary  condition  is  the  consequence  of 
the  fact  that  the  flux  at  the  interface  has  to  be  equal  to 
the  rate  of  heterogeneous  electrode  reaction,  i.e.: 

xm0  *  -  hCcV-  (32) 

where  superscript  "o"  denotes  the  surface  condition.  To 
satisfy  the  stoichiometric  relation  of  equation  (23),  one 
has,  at  x=0 : 


-  - 

f  \ 

_  ^2 

3x 

1 

o 

II 

X 

^2  -v 

3x 

J 

This  is  the  second  boundary  condition.  The  ct.-.er  two 
boundary  conditions  are: 
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3  •  at  x=oo  ,  C^  =  C° , 

•u 

4.  at  x=co  ,  C2=  C2  • 

That  is,  the  concentrations  far  away  from  the  electrode/ 
electrolyte  interface  do  not  change. 


The  flux  in  equation  (32)  is  proportional  to  the  curren 
density  as  a  function  of  time: 
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CHAPTER  4 

EXPERIMENTAL  WORKS 

A  potential-step  was  generated  from  a  function 
generator  and  was  applied  to  the  working  electrode/elec- 
trolyte  interface  in  the  electrochemical  cell.  The  ca¬ 
thodic  current  was  thus  sampled  and  stored  in  memory  at 
a  constant  time  interval  by  a  microcomputer.  The  resulting 
current/time  transient  was  then  plotted  by  a  x-y  recor¬ 
der  or  printed  out  by  a  line  printer  by  recalling  these 
current  data  from  the  memory  of  the  microcomputer. 

1  The  Electrochemical  Cell 


The  current  transient  under  a  potential  step  con¬ 
dition  was  obtained  in  a  three-electrode  cell.  The  cell 
and  electrodes  configuration  is  sho-wn  in  Figure  1.  The 
working  electrode  was  a  nickel  microelectrode.  It  was  made 

of  a  nickel  wire  which  was  pressure-fitted  into  a  teflon 

.  2 
cylinder.  The  flat,  exposed  surface  was  0.013  cm'" .  A 

saturated  calomel  electrode  (SCE)  and  a  cadmium  bar  were 
used  as  the  reference  and  che  counter  electrodes,  res¬ 
pectively  . 


The  cell  was  filled  with  electrolyte  of  either 
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cadmium  nitrate  solution  or  potassium  nitrate  solution. 

The  pH  value  of  the  solution  was  controlled.  Potassium 
chloride  was  used  as  the  supporting  electrolyte.  The  electro¬ 
lytic  solution  was  not  stirred  during  the  experiment,  so  that 
the  conditions  of  semi-infinite  linear  diffusion  were  maintained. 

Two  experimental  parameters  were  varied: 

1)  the  ion  concentration  of  the  solutions  and 

2)  the  magnitude  of  the  potential-step  applied  to  the 
interface. 

4.2  The  Instruments  and  Electrical  Setup 

The- block  diagram  for  the  experiment  setup  is  shewn  in 
Figure  2.  The  instruments  shewn  in  Figure  2  contain  the 
following  equipment: 

a)  A  Princeton  Applied  Research  (PAR)  Model  1'~’3  Potentiostat/ 
Galvan os tat :  In  the  experiment,  the  operating  mode  was  set 
at  CONTROL  E. ,  which  means  that  the  instrument  functioned 
as  a  potentiostat.  The  current  was  measured  -while  keeping 
the  potential  of  the  working  electrode  constant.  This  was 
accomplished  by  setting  the  counter  electrode  to  the 
required  level  so  that  the  'working  electrode  potential  was 
at  a  programmed  potential  with  respect  to  the  reference 
electrode.  The  instrument  is  provided  with  a  Me  del  1"6 
Current- to -Voltage  Converter  which  is  the  "basic"  plug-in 
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module.  The  R1  176  provides  a  do  voltage  output  which  is 
proportional  to  the  current.  The  cell  current  can  be  eitner 
displayed- on  the  panel,  by  a  x-y  recorder,  or  monitored  by 
the  microcomputer.  Connection  to  the  cell  was  made  through 
the  external  cable.  The  counter  electrode  was  connected  to 
the  red  clip,  the  working  electrode  to  the  green  clip  and 
the  reference  electrode  to  the  Electrometer  Probe  which  has 
a  very  high  impedance,  thereby  insuring  us  that  the  current 
was  flowing  to  the  reference  electrode, 

b)  A  PAR  Model  175  Function  Generator!  This  is  a  programmable 
waveform  generator  which  has  two  operating  modes,  SWEEP  and 
PULSE.  The  former  generates' a  sequence  of  triangular  waves, 
while  the  latter  generates  a  sequence  of  step  function,  Eor 
the  potential  step  experiment,  the  operating  mode  was  set- a 
PULSE  mode  and  the  pulse  width  selector  at  STEP  position. 
The  magnitude  of  the  step* function  was  set  by  the  setting  o 
the  3  potential  (upper  limit)  in  the  POTENTIAL  section  on 
the  front  panel  of  this  instrument. 

c)  An  Intel  CPU  S080A  Bases  Microcomputer  ( z^K  PA.',.):  The  micro 
computer  was  the  central  part  of  the  experimental  setup. 

The  functions  of  the  microcomputer  in  the  experiment  were: 

(1)  to  trigger  the  Function-* Generator  which  activated  the 
potential  step  presented  earlier, 

(2)  ...to  sample  the  current  outputs  from  the  M  17c  Converter 

- and  store -them  in- the- memory .  — 

(3)  to  convert  these  digital  data  back  into  analog  signals 
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and  then  plot  them  on  the  x-y  recorder. 

(4)  to  print  the  digital  data  on  the  line  printer. 

(5)  to  store  these  digital  data  in  the  tapes. 

The  peripheral  equipment  of  the  microcomputer  include  a 
CRT  (TV  screen),  the  keyboard,  a  cassette  recorder  and  a 
line  printer. 

d)  A  Houston  Model  RE  0074  x-y  Recorder:  The  Recorder  receives 
the  analog  outputs  from  the  DAC  (Digital-Analog  Converter) 
and  plots  them  on  the  graph  paper. 

4.3  The  Computer  Program 

A  computer  program  ’written  in  BASIC  language  -was  used  no 
carry  out  the  experiment.  The  program  consists  of  a  main  program 
and  a  machine-language  subroutine  called  "MUG".  This  subroutine 
includes  an  execution  statement  which  generates  a  trigger  signal 
to  activate  the  Function  Generator.  A  potential  step  was  then 
immediately  applied  to  the  working  electrode,  and  thereafter 
the  computer  started  to  sample  and  store  the  current  output  at 
a  fixed  time  interval.  The  analog  signal  was  converted  into  a 
digital  datum  by  a  12  bits  ADC  (Analog-Digital  Converter)  and 
then  stored  in  the  memory.  When  the  number  of  samples  reached 
a  preset  value,  the  sampling  routine  was  terminated.  The  digital 
signals  stored  in  the  memory  were  binary  numbers.  These  data 
were  converted  into  decimal  numbers  and  then  converted  into 
analog  signals,,  and  finally  plotted  on  the  x-y  recorder  and 
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printed  on  the  line  printer.  This  was  executed  by  the  main 
program.  The  main  program  and  the  "MUG"  subroutine  are 
presented  fully  in  Appendix  A1  and  A2 . 

4.4  Experimental  Conditions  and  Procedures 


In  the  experiment,  the  initial  pH  value  of  the  electro¬ 
lyte  was  maintained  at  J.O,  which  was  the  condition  used  in  the 
actual  electrochemical  impregnation  process.  This  was 
accomplished  by  titrating  the  electrolytic  solution  with  a 
concentrated  HC1  solution.  The  concentrations  of  nitrate 
ion  used  were  0.005M,  0.002574  and  0.00125M  prepared  from  either 
the  KNC^  or  Cd(NO^)  solution.  The  cadmium  ion  concentrations 
were  0.0025M,  0.00125M  and  0.00Q625M,  All  the  electrolytic 
solutions  contained  0.5M  KCi  as  the  supporting  electrolyte. 

It  was  .found  later  that  the  electrolysis  of  water  made 
significant  contribution  to  the  total  current,  because  the 
current  obtained  was  significantly  higher  than  the  limiting 
current.  In  order  to  obtain  the  current  due  only  to  the  reduction 
of  nitrate  ions,  several  .runs  with  solutions  which  contained  no  . 
nitrate  ion  were  conducted.  The  solutions  used  for  this  purpose 
were  either  water  or  C.0G25M  CdCl^  solution.  Beth  solutions 
contained  C.5M  KCi  as  the  supporting  electrolyte.  This  current 
was  subtracted  from  the  current  obtained  in  t he  presence  of 
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nitrate  ions  in  the  solution.  The  resulting  current  is  the 
faradic  current  of  the  reduction  of  nitrate  ions. 

The  magnitudes  of  the  potential-step  used  in  this 
study  were  -0740V,  -0 . 60V  and  -0 . 80V  versus  the  equilibrium 
potential  which  was  at  approximately  -0.35^7  vs  SCE, 

The  experimental  procedures  are  listed  blows 

(1)  The  electrical  circuit  was  set  up  as  shown  in  Figure  2, 
except  that  the  cell  was  disconnected  from  the  Potent iostat/ 
Galvanostat  by  setting  the  cell  selector  on  the  instrument 
to  the  OFF  position.  A  proper  current  range  was  then 
chosen  from  the  Current  Range  Switch  to  make  sure  that 

the  current  output  would  not  be  overloaded.  The  Function 
Generator  was  initialized  by  depressing  the  INITIAL 
control  pushbuttom.  The  initial  potential  was  set  at 
zero  volts. 

(2)  A  known  volume  of  electrolytic  solution  and  an  equal 
volume  of  supporting  electrolyte,  KC1,  were  added  to  the 
cell  and  titrated  with  HC1  solution  to  a  pH  value  of  3*0. 

The  solution  was  then  deaerated  by  bubbling  purified 
nitrogen  through  the  stirred  solution  for  about  ten 
minutes.  The  gas  continued  to  pass  above  the  solution 
during  one  experiment. 

(3)  After  the  working  electrode  surface  was  polisr.ed  by  a 

3/0  alumina  paper,  it  was  rinsed  thoroughly  with  distilled 
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and  deionized  water,  and  then  was  positioned  in  the  cell. 

The  electrodes  (the  working  electrode,  counter  electrode 
and  reference  electrode)  were  then  connected  to  the 
Potentiostat/Galvanostat  as  described  in  section  4. 3 • 

After  the  electrodes  were  correctly  connected,  the  cell 
selector  was  switched  to  the  EXTERNAL  CELL  position.  The 
rest  potential  of  the  solution  could  be  determined  by 

turning  the  knob  in  the  APPLIED  POTENTIAL/ CURRENT  section 
of  the  front-panel  of  the  Potentiostat/Galvanostat  so 
that  the  current  reading  displyed  on  the  front  panel 
meter  was  zero.  The  applied  potential  was  then  the  rest 
potential. 

(4)  The  upper  limit  (B  potential)  on  the  Function  Generator 
was  set.  The  value  of  the  -3  potential  plus  the  rest 
potential  was  the  total  potential  applied  to  the  electror 
chemical  cell. 

(5)  The  experiment  was  initiated  by  running  the  computer 
program.  The'  current/time  data  were  stored  in  the  computer’s 
memory  and  -then  the  results  were  converted  to  line 

printer  output. 

4. 5  Experimental  Results 

Two  sets,  a  total  of  twenty-two  runs,  were  made.  The 
first  set  of'  runs  was  conducted  in  solutions  which  contained 
no  cadmium  ion.  The  second  set  of  runs  was  made  in  solutions 


34 


which  contained  cadmium  ions.  For  each  set  of  runs,  clank  runs 
were  made  in  the  solutions  without  the  presence  of  nitrate 
ion  for  the  purpose  of  determining  the  background  current 
aue'to  the  water  decomposition  reaction.  The  solution  for 
blank  runs  for-the  first  set  was. water,  while  the  solution 
for  blank  runs  for  the  second  set  was  made  by  using  the 
0.0025M  CdC^  solution. 


The  current-density-vs .-time  data  printed  by  the  H  14 
Line  Printer  are  shown  in  Table  3.1  through  Table  B.22  in 
appendix  3.  The  measured  current  was  converted  to  the  current 
density  which  appears  in  those  tables.  This  is  calculated 
by  dividing  the  current  by  the  electrode  surface  area 
0.013ca2. 


Tables  3.1  through  3.9  are  the  first  data  set  which  used 
potassium  nitrate  as  the  electrolytic  solution,  while  Tables 
3.10  through  3,12  are  the  background  currents  of  water  used 
to  subtract  from  the  first  data  set. 


Tables  3.13  through  3.21  are  the  second  data  set  which 
used  cadmium  nitrate  as  the  electrolytic  solution,  while 
Table  3.22' is  the  background  current  of  water  in  one  presence 
of  cadmium  ions. 


A  series  of  current-density-vs , -time  curves  witn  various 
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electrolytic  solution  concentrations  and  applied  overpotentiais 
are  shown  in  Figures  3  through  6.  In  these  figures,  the 
current  due  to  the  electrolysis  of  water  has  been  eliminited 
by  using  the  routine  tabulated  in  Table  1. 

Figures  3  through  5  are  the  current-density-vs . -time 
curves  obtained  in  the  solutions  which  did  not  contain  cadmium 
ion.  Figure  3  shows  the  current-density-vs . -time  curves 
obtained  in  a  solution  containing  0.005M  nitrate  ion  and  at 
potentials  -0.40V,  -O.oOV  and  -0.80V  from  the  equilibrium 
potential.  Figure  4  shows  the  current-der.sity-vs , -time  curves 
obtained  in  a  solution  containing  0.0025M  nitrate  ion,  and  at 
potentials  -0.40V,  -0.60V  and  -0.807  from  the  equilibrium 
potential.  Figure  5  shows  the  current-density-vs . -time  curves 
obtained  in  a  solution  containing  0.00125M  nitrate  ion  ax  -0.40V, 
-0.60V  and  -0„80V_from  the  equilibrium  potential.  Due  to  the 
lack  of  blank  runs  with  0.00125M  and  0.000625M  cadmium  ion  con¬ 
centrations  and  -0.60V  and  -0.80V  from  the  equilibrium  potential, 
only  one  current-density-vs . -time  curve  was  obtained  with  0.0G5M 
nitrate  ion  and  0.0025M  cadmium  ion  concentraxion  at  the  poten¬ 
tial  of  -0.40V  from  the  equilibrium  potential.  This  is  shown 
in  Figure  6. 

6  Discussion  of  Zxxeriner.xal  Results 


Some  characteristics  can  be  seen  from  the  curves  shown 
in  Figures  3  through  6. 
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Table  Is  Evaluation  of  Currer.t-Density-vs .  -Time  lata 
Plotted  in  Figures  3  Through  6. 


Figure  is 
plotted 

by  subtracting  the 
corresponding  current 
density  at  each  time 

Doint  in 

from  that 

of  in 

(1) 

. 

Table  B.l 

Table  3.10 

Figure  3 

(2) 

Table  3.2 

Table  3.- 11 

(3) 

Table  3.3 

Table  3. 12 

(1) 

Table  3.4 

Table  3.10 

Figure  4 

(2) 

Table  3.5 

Table  3.11 

(3) 

Table  3.6 

Table  B.12 

(1) 

Table  3.7 

Table  3.10 

Figure  5 

(2) 

Table  3.8 

Table  3.11 

(3) 

Table  3.9 

Table  3.12 

Firstly,  the  current  density  decays  with  time  in  a  ] 

hyperbolic  fashion.  This  is  the  result  of  depletion  of  the  j 

iN'O^  ions  near  the  electrode  surface.  Theoretically,  there  j 

is  an  initial  sharp  rise  of  current  associated  with  double-layer 
(27) 

charging  ,  due  to  the  application  of  a  potential-step .  This  was  \ 

not  observed  in  the  experiment,  however,  because  the  time  for  j 

charging  is  usually  very  short  as  compared  to  the  electrolysis  1 

time. 

;i 

i! 

Secondly,  the  sharp  current  drops  in  the  beginning  indicate  j 

that  very  sharp  concentration  gradients  are  established  for  the  j 

:  I 

nitrate  ions  and  the  hydroxy  ions  as  soon  as  the  surface 

reaction  (23)  takes  place.  The  concentration  gradients  are  ! 

i 

gradually  decreased  due  to  the  diffusion  layer  extending  in  ^ 

the  direction  of  decreasing  concentration  gradients.  The  i 

result  is  that  the  decay  of  the  current  density  is  not  as  fast  j 

as  in  the  beginning. 


The  curves  for  the  case  of  potential  at  -0.80V  from  the 
equilibrium  potential  and  solutions  containing  0.G025M  and 
0.00125M  nitrate  ion,  as  shown  in  the  upper-most  curves  in  Figures  - 
and  5»  respectively,  are  somewhat  different  from  t ns  others. 

These  two  curves  do  not  follow  the  same  trends  as  their  lower 
potential  counterparts.  It  is  suggested  that  the  reactant, 
nitrate  ion,  near  the  electrode  surface,  is  used  up  in  a  very 
short  time  and  another  reaction  must  be  taxing  place.  This 


behavior  is  the  mos*  pronounced  as  shown  in  the  upper-most 
curve  in  Figure  5*  A  maximum  current  is  observed  at  time 
equal  to  about  0.5  msec.  This  curve  thus  strongly  supports 
our  explanation  that  some  other  reaction  is  taking  place  at 
that  time.  On  the  other  hand,  the  absence  of  this  behavior 
in  the  curve  for  the  case  of  O.OOyM  nitrate  ion  shown  in  the 
upper-most  curve  in  Figure  3  suggests  that  the  nitrate  ion 
concentration  near  the  electrode  surface  is  not  zero. 

Two  parameters  were  varied  in  the  experiment,  i.e., 
the  concentration  of  the  electrolyte  solutions  and  the  value 
of  the  potential  from  the  equilibrium  potential  (rest  potential). 
The  changing  of  either  the  bulk  concentration  of  the  electrolyte 
solutions  or  the  overpotential  that  is  applied  to  the  working 
electrode  will  change  the  value  of  the  current  density.  For 
a  given  bulk  concentration,  the  increase  in  the  applied  over- 
potential  increases  the  electrochemical  reaction  rate  on  the 
electrode  surface.  The  effect  of  the  applied  overpotential  on 
the  current  density  can  also  be  observed  in  Figures  3  through  5> 
In  each  of  these  figures,  a  plot  of  the  limiting  current  - 
density- vs .-time  curve  is  also  shown.  The.  limiting  current 


density 

is  calculated 

by  equation  (9) 

the  nitrate  ion. 

•with  the 

corresponding 

bulk  concentration  of 

Fecal I 
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current 

is  defined  as 

maximum  current 
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overpoter.tial 

is  increased.  T 

he  surfac 

e  concentration. 

of  the  reactant  for  the  case  of  limiting  current  is,  in  fact, 


^3 


zero  at  any  tine.  After  examining  Figures  3  through  5 >  we  find 
that,  for  the  case  of  -0.30V,  some  of  the  values  of  current 
density  exceed  the  limiting  current  density.  This  phenomenon 
becomes  more  significant  when  the  bulk  concentration  of  pota¬ 
ssium  nitrate  is  decreased.  The  reason  for  this  phenomenon 
probably  is  that  at  the  high  applied  overpotential ,  many 
reactions  can  subsequently  be  expected  to  take  place.  The 
maximum  current  appearing  in  the  curve  for  the  case  of  -0.50V 
is  a  result  of  this  behavior. 

When  the  applied  overpotential  is  fixed,  the  hetero¬ 
geneous  reaction  rate  constants  in  both  the  cathodic  and  anodic 
directions,  K,  and  ,  are  fixed.  The  supply  of  the  nitrate 
ions  on  the  electrode  surface  comes  only  in  the  way  of  diffusion 
of  the  nitrate  ions  from  the  bulk  solution.  Increasing  the 
bulk  concentration  will  increase  the  rate  of  nitrate  ion  supply 
to  the  electrode  surface,  thusly  increasing  the  current  density. 
Figures  7  through  9  show  the  effect  of  the  bulk  concentration 
of  the  electrolytic  solution  on  the  current  density  with  fixed 
applied  overpotential.  From  the  three  figures,  we  can  see  that 
the  current  density  is  approximately  proportional  to  the  bulk 
concentration  of  potassium  nitrate,  except  for  the  case  of  -C: . 

When  electrolytic  solution  contains  cadmium  ions,  a 
different  set  of  results  was  obtained.  The  experimental 
current-density-vs . -time  curve  obtained  in  a  0.00251-! 
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Effect  of  the  Bulk  Concentration  of  the  Electrolytic  Solution 
on  the  Current  Density  with  Fixed  Applied  Overpotential . 
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Cd(N0^)2  solution  and  at  -0.40V  applied  overpotential  is  shown 
in  Figure  6.  In  addition  to  the  electrochemical  reduction  of 
the  nitrate  ions  on  the  electrode  surface,  a  precipitation 
reaction  also  takes  place  in  the  solution  because  the  cadmium 
ions  in  the  electrolytic  solution  can  coprecipitate  with  the 
hydroxy  ions  produced  by  the  electrochemical  reaction  on  the 
electrode  surface.  The  current  density  is  expected  to  be 
greater  than  that  when  the  cadmium  ion  is  absent  from  the 
electrolytic  solution.  Figure  10  shows  the  current-density - 
vs. -time  curves  obtained  in  a  0.005M  KNO^  solution  and  in  a 
0.0025M  Cd(N0^)2  solution.  The  applied  overpotential  is  -0.40V 
for  both  cases.  From  the  curves,  one  can  see  that  the  current 
density  obtained  in  the  0.0025M  Cd(N0^)2  solution  is  larger 
than  that  obtained  in  a  0.005M  KNC^  solution.  The  reason  for 
this  is  that  the  consumption  of  the  hydroxy  ions  due  to  the 
precipitation  reaction  establishes  a  driving  force  to  produce 
more  hydroxy  ions  and  thus  promote  the  electrochemical  reaction. 
Some  interesting  facts  may  be  seen  from  the  figure.  Firstly, 
the  difference  between  two  current  densities  increases  as  time 
increases,  finally  reaching  a  constant  difference;  secondly, 
the  rate  of  decrease  of  the  current  density  with  time  is  slower 
for  the  case  of  using  Cd(N0^)2  as  the  electrolytic  solution. 

One  may  thus  obtain  a  larger  steady  state  current  density  using 
Cd(N0^ ) 7  as  the  electrolytic  solution. 
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CHAPTER  5 

NUMERICAL  SOLUTIONS 
OF  THE  THEORETICAL  MODEL 

5 • 1  Orthogonal  Collocation  Method 

There  are  many  numerical  techniques  which  can  be  used 
to  solve  the  differential  equations  of  the  theoretical  model 
presented  earlier  in  Chapter  3* 

Orthogonal  collocation  method  was  selected  as  the  method 
of  solution.  The  orthogonal  collocation  method  requires  very 
little  implicit  mathematics  and  results  in  a  large  savings  in 
computer  time  over  the  common  finite  difference  scheme  such  as 
the  Crank-Nicholson  method.  A  brief  description  of  the  ortho¬ 
gonal  coxlocation  method  is  presented  in  Appendix  C. 

5*2  Working  Equations 


5*2.1  Dimensionless  Form  of  the  Differential  Equations 


The  orthogonal  collocation  method  requires  the  region  of 
solution  to  be  restricted  to  the  interval  C C , 1 3  .  This  is 
accomplished  by  normalizing  the  spatial  variable  by  the 
parameter ,  £ ,  the  diffusion  thicxr.ess,  ^  =x/ /  . 


dimensionless  forms s 

* 

p  —  ~  •  r*  _  - 

U  ,  t  L>  ~  — 

1  Cl  2  c? 

1  1 

the  dimensionless  time  is  defined  as: 


Substituting  the  above  dimensionless  variables  into  equations 
(17)  and  (19)  and  rearranging  leads  to  the  following  dimension¬ 
less  eauations: 


.(.Zulh 
1  vay 


*  2„  * 

^  ~,2  *  „  *q  ZK-v 

•2T  ‘O'?2  ■  r  1 


where  K=  — 


and  DKsp=  -^52- 


are  the  dimensionless  homogeneous  reaction  rate  constant  and 
solubility  product,  respectively.  Similarly,  equations  (21) 
and  (22)  become: 


_  _2„  * 
U1  \  1>  U1 


*  ?  * 
_  ycz 


V>J 
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The  initial  conditions  are: 


1.  at  r  =0 ,  CL  =1, 


for  all  0 


O  s  l- 


2.  atr« 


=0,  C2~=  C°/C°  ,  for  all  0  <  y  < 


1. 


The  first  two  boundary  conditions  come  from  the  dimensionless 
form  of  equations  (32)  and  (33)* 


D1  C°  f  >ci 

1  •  (  -l—i-  )  I  0  1 


*  N 


h  fay  J$=0 


k2  (c'J)3  <C°V  -  (X.CJ)  c°’ 


=  0  , 


3i  f3c-  1 

2.  3(  — )  %-£- 

C2  L23 


f  3C; 


=  0 


:=o 


The  other  two  boundary  conditions  are: 


3 .  at  ^  =1 ,  C1  =1  , 


for  X  >  0  , 

4.  at^-**  C2*=C^/cr,  for  *£•  >  0 . 


5.2.2  Discretized  Equations 


The  solution  is  approximated  by  a  (L’+2 )  th-order  Legendre 
polynomial  which  satisfies  the  differential  equations  and 
boundary  conditions  exactly  at  N+2  points.  Those  points  are 
chosen  to  be  zeros  of  the  Legendre  polynomial  of  degree  over 
the  interval  (0,1)  and  with  the  boundary  points  0  and  1. 
Discretizing  the  spatial  derivatives  at  those  points  leads  to 
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the  following  set  of  first-order  ordinary  differential  equations 


i=  2, . N-rl,  (39) 


i  3  2 , .  N + 1 ,  ( IrO ) 

where  ^  ^ ,  i=2,...  N+l  are  the  N  inferior  collocation  points 
while  1  and  N+2  are  exterior  collocation  points  which 
correspond  to  the  boundary  points  at  ^  =0  and  =1, 
respectively,  and  0=  ^  %  ^  ^  N+2  =1,  i=2,  .  . .  N  +  l.  3j_  , 

is  a  (N+2)  by  (N+2)  coefficient  matrix  which  cones  from  the 
discretization  of  the  second  derivative  (  d2ci*/d^  or 
O  C2  / S  *5  )  eacn  collocation  point.  A  detailed  explanation 
of  matrix  3  is  shorn  in  Appendix  C.  C,  .,'£■)  and  C-,  (t  .,t) 
are  the  concentrations  of  species  1  and  2,  respectively,  at 
fine  'C  and  position  ^  ^ . 
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where  ^  is  an  element  of  another  (N+2)  by  (N+2)  coefficient 
matrix  which  comes  from  the  discretization  of  the  first-order 
derivative  (  2>C^  or  "^C2  )  at  each  collocation  points . 

Since  only  boundary  conditions  1  and  2  (equations  (43)  and  (44)) 
require  the  evaluation  of  the  first-order  derivative  at  *?  1 =0 , 
only  the  first  row  of  the  matrix  A,  i.e.,  A,  is  used  here. 

1  I  o 

The  initial  conditions  for  both  cases  are  the  same: 

I.C.  ’s 

1.  at  -f=0,  c1*(^i<T)=  1>  i  =  l....N+2, 

2.  at  ^  =0,  )=  C2//Ci:  ’  i=lf.M+2. 


5*2.3  Calculation  of  Current  Density 


The  current  density  is  related  to  the  flux  of  the  nitrate 
ion  at  the  electrode  surface  by  equation  (34).  The  discretized 
and  dimensionless  form  of  equation  (34)  is; 
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3 . 3  The  Introduction  of  a  Saline  Point  to  the  Discretized 
Equations 


The  diffusion  thickness  introduced  in  the  derivation 
of  the  dimensionless  form  in  the  previous  section  can  he 
justified  by  selecting  the  total  electrolysis  time  such 
that  diffusion  effects  are  negligible  at  x=£  during  that 
period.  At  the  beginning  of  the  electrolysis,  the  diffusion 
thickness  selected  according  to  the  total  electrolysis  time  is 
too  large.  This  leads  to  a  concentration  profile  which  is 
flat  in  most  of  the  interval  with  a  sharp  drop  in  a  very  small 
region  near  *5-  =0.  This  would  require  a  large  number  of 
collocation  points  leading  to  a  very  stiff  set  of  ordinary 
differential  equations. 

This  difficulty  can  be  overcome  by  the  use  of  a  method 

f  2  Q  \ 

called  spline  collocation'  '  based  on  the  diffusion  boundary 
concept.  This  method  maintains  a  fixed  low  number  of  collo¬ 
cation  points  in  the  interval  which  is  very  small  initially 
and  will  be  increased  as  time  increases  to  account  for  the 
thickening  of  the  diffusion  layer.  In  the  regions  outside 
this  interval,  the  concentration  is  fiat.  As  a  natter  of 

fact,  the  spline  point  can  be  located  as  close  to  the  inter- 

(  ?  q  \ 

face  as  to  achieve  any  desired  accuracy  , 

t  the  concentration  gradient 


In  order  to  be  sure  tha 


at  the  spline  point  is  zero,  the  concentration  at  the  Hth 
collocation  point,  i.e.,  the  last  one  before  the  spline 
point,  is  tested  to  assure  that  it  is  within  a  very  small 
range  of  the  bulk  condition.  At  the  time  when  such  a  test 
fails,  the  spline  point  is  moved  further  into  the  solution, 

Introducing  a  spline  point,  ^  s ,  G<  ^  g  4.  1 ,  a  new  variable, 
z= 'V /  ^s,  is  introduced.  The  discretized  equations  we  obtained 
for  the  case  of  solution  containing  cadmium  ions  are: 


dC 

1 

df  It; 


— —  r  :;+2  *  -l 

*  <  v2  >  <  r*  )  n  3<,-h  <-”T> 

Lfe  “  •  J 


1=2 , «  e  9  N  + 1 1 


i 

r 

■  (  >2  ) 

1  y  3,  ,c  (2,,r) 

-2K  f 

2.  S3 

^  *  ■  ■  -  f  j  j  j 

j=i 

l 

^C2  (zL,X)) 


o  "0 

*  2 


1=2 , . . .n*i . 


For  the  case  of  solution  containing  no  cadmium  ion, 


c  n 

-/  i 


Equations  (43)  through  (46),  which  cone  from  the  four  c cur. car; 
conditions,  become: 

D  cb  ?I+2 

la  (“fs  }  (  Al,jCl  (sj»/r[)  +  |^2'C1)3c2 

-  (X-)c;(t1>T)]  =0,  for  allt>0.  (32) 


2.  3(- 


Di  .  <  ?il_2 


L>[Z  Wi'fyX)]  ♦  [E  ai,;c2*^.t)] 

2  u  j  =  l  ^  ^  -  =  l  > 


for  all  X>°* 


3-  ci*(^s.t)“  ci*(::;.2'tl=  1.  ^or  ail  T>0, 

4-  52*^S'T>*  s2#‘*a»2»t)*«I/*i.  --i  Xfi. 


The  initial  conditions  become: 
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A 


4 


T  r  '  <3  • 

1  •  •  o  * 

* 

1.  at  ^=0,  C1  (zi,T)=  1.  i=l,...N+2, 

2.  at  r=0,  C2*(zi>r)=  C2/C2,  i=l,...N*2. 


The  current  density  should  he  changed  to  s 


5.4  Solution  Procedures 


For  the  case  of  a  solution  containing  cadmium  ions, 
equations  (48),  (49)  and  (52)  through  ( 55 ) »  with  the  initial 
conditions,  are  the  working  equations  used  to  solve  for  the 
concentration  profiles  of  N 0 ^ "  and  OH-  ions.  The  Cd"1"^  ion 
concentration  is  calculated  by  the  electroneutrality  conditic 
in  the  solution: 
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For  the  case  of  a  solution  containing  no  cadmium  ion, 
equations  (50),  (51)  and  (52)  through  (55)  with  the  initial 
conditions,  are  the  working  equations. 

Equation  (56)  is  used  to  calculate  the  current  density 
at  each  time  /f  for  both  cases. 

5 • 5  Computer  Program  Structure 


The  computer  program  for  solving  the  discretized 
equations  consists  of  a  main  program  and  5  main  subroutines. 
The  first  main  subroutine  evaluates  the  collocation  points 
(roots)  of  the  Jacobi  polynomial  of  order  N,  as  well  as  the 
firs  and  second  derivatives  tne  polynomial  at  the  roots. 
These  derivatives  are  then  used  in  a  second  main  subroutine 
to  caxculate  the  discretization  coefficient  matrices  A  and  3. 
The  third  main  subroutine  is  used  to  perform  a  semi-implicit 


integration  using  Gear's  routine 


(29) 


Inis  subroutine  cans 


4  external  subroutines.  The  first  of  these  contains  the 

explicit  expressions  (48)  and  (49)  (or  (50)  and  (51)  in 

another  case)  for  the  discretized  coupled  first  order  differentia- 

equations  containing  the  3.  .  terms;  the  second  contains  the 

-  >  j 

Jacobian  matrix  for  the  same  equation;  the  third  performs  the 
first  stage  of  Gaussian  elimination^  ;  of  the  Jacobian  matrix; 
tr.e  _ast  performs  the  back  substitut-on  of  Gaussian  elimination 


of 
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of  the  Jacobian  matrix.  The  first  external  subroutine  called 
by  the  integration  subroutine  also  calls  another  subroutine 
coming  from  the  IMSL  package.  This  subroutine  solves  equations 

(52)  and  (53)  simultaneously  to  give  the  surface  concentrations 

*  * 

of  C.  and  C-  .  The  IMSL  subroutine  calls  another  function 

1  c. 

which  feeds  the  equations  to  be  solved.  Before  writing  the 

*  * 

function,  (surface  concentration  of  Cp  )  is  solved 

-g.  -fr 

in  terms  of  (z-,^),  i=l ,  .  .  .'i-^2 ,  and  (z-.'T),  i  =  l ,  .  .  ..'iT2  , 

using  equation  (53) »  and  then  substitute  into  equation  (52) 

* 

to  eliminate  (z^,^).  Thus,  only  one  equation  which 

* 

contains  only  one  unknown,  is  needed.  Cnee  we 

obtain  (z^.'-j-),  C?  (z^,^)  can  be  calculated  easily. 

The  Jacobian  matrices  for  both  cases  can  also  be  evalu- 

* 

ated  after  substituting  the  expression  of  (s-,,-^)  into 
equations  (43)  and  (49)  (or  (50)  and  (51)  in  another  case). 

For  the  case  of  solution  containing  no  cadmium  ion,  the  Jacobian 
matrix  is: 
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which  is  a  2N  by  2N  matrix.  For  the  case  of  solution 
cadmium  ions,  the  Jacobian  matrix  is  the  same  except 
the  terms  in  the  diagonal  line  of  the  third  quadrant 
subtracted  by  X  and  the  terms  in  the  diagonal  line  of 
fourth  quadrant  are  all  subracted  by  2K( f+( TDXsp/C^  v 


containing 

that 


are  all 
the 


T) 


The  values  obtained  from  the  integration  subroutine  are 
then  used  in  a  fourth  main  subroutine  to  calculate  the  current 
density  by  equation  (56),  The  last  main  subroutine  is  used  to 
evaluate  the  concentration  at  any  desired  point  between  0  and  1 
by  interpolation  using  the  values  at  the  collocation  points. 

The  structure  of  the  whole  program  is  shown  in  Figure  11. 
Programs  from  both  cases  are  listed  in  Appendix  D.  Programs 
herein  were  executed  on  an  AMDAHL  4?0/V6  computer. 
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CHAPTER  6 

FITTING  OF  THE  THEORETICAL 
MODEL  WITH  EXPERIMENTAL  DATA 

6 . 1  The  Determination  of  the  Electrochemical  Reaction 
Kinetics  Parameters 

The  unknown  heterogeneous  reaction  rate  cons¬ 
tants  which  appear  in  equation  (52)  were  determined 
from  the  current-der.sity-vs ,  -time  data  obtained  in  a 
solution  which  did  not  contain  cadmium  ion.  The  di¬ 
fferential  equations  (50)  and  (51)  which  describe  the 
transport  process  were  solved  by  computer  to  determine 
the  rate  constants.  Comparing  equation  (31)  with  equa 
tion  (31a),  one  can  see  that  both  K1  and  K9  depend  on 
the  applied  potential  with  exponential  dependence. 
Thus,  X.  and  ^  should  be  constant  when  the  applied 
potential  is  maintained  at  a  constant  value  during  the 
electrolysis.  K1  and  K^  were  obtained  by  comparing  the 
experimental  data  in  the  absence  of  cadmium  ion  in  the 
electrolyte  with  the  computer-calculated  current  der.si 

The  experimental  current-der.sity-vs  , -time  data 


data 


with  the  electrolyte  containing  0.0057  nitrate  ions 
and  the  applied  overpotential  of  -0,407  was  first 
fitted  to  the  computer-simulated  data  to  obtain  the 
K1  and  values.  Due  to  the  lack  of  information 
about  the  values  of  X1  and  ,  rough  estimates  of  X. 
and  K0  were  evaluated  by  substituting  the  applied 
potential  V,  which  is  equal  to  the  value  of  the  applied 
overpotential  plus  the  rest  potential,  to  equations 
(31)  and  (31a).  The  estimated  values  of  X,  and  X,  were 
increased  or  decreased  until  the  computer-calculated 
current-density-vs .-time  curve  had  the  same  shape  as 
the  same  curve  obtained  experimentally,  (note,  the 
background  water  decomposition  current  has  been  sub¬ 
tracted.  ) 

After  we  obtained  the  right  shape,  the  anodic 
reaction  rate  constant,  X^ ,  was  then  fixed.  By  changing 
X, ,  one  can  make  the  entire  curve  shift  upward.  Con¬ 
versely,  decreasing  the  X,  value  will  decrease  the 
current  density  at  each  time  point  and  thus  shift  the 
entire  curve  downward.  A  set  of  current-density-vs . - 
time  curves  with  fixed  X9  values  and  different  values 
of  X.  are  plotted  in  Figure  12.  The  values  of  X,  and 
X,  when  the  concentration  of  nitrate  ions  is  0.0057, 
and  the  applied  overpotential  is  -C.mJV  w as  obtained 


Electrolysis  Time,  ms 
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by  interpolating  the  experimental  curve  b etween  chose 
curves  in  Figure  12.  The  values  of  X1  and  X 2  which 
achieved  the  best  approximation  of  the  experimental 
data  are  K^=  0.12  X  10"^  and  X2  =  0.6  xlO--5,  respectively 

These  values  were  then  used  to  predict  the  cu- 
rrent-density-vs .-time  curves  at  the  other  ion  con¬ 
centrations.  Figure  13  shows  the  computer-simulated 
current  density  curves  and  experimental  curves  at  va¬ 
rious  nitrate  ion  concentrations . 

The  values  of  X^  and  X2  at  the  applied  overpo¬ 
tential  of  -0.60V  were  different  from  those  at  -0.40V. 
The  same  procedure  was  used  in  evaluating  the  values  of 
and  X2  at  the  condition  of  -O.cGV.  Figure  14  shows  a 
set  of  calculated  current-density-vs . -time  curves  with 
fixed  K2  valued  and  different  values  of  X., .  The  value  o 
X*  at  the  condition  of  -0,60V  was  then  obtained  in  the 

i. 

same  way  as  before.  The  values  of  X,  and  X0  were  deter- 
mined  to  be  0.10  x  10  J  and  0.50*10  ,  respectively. 

These  values  were  then  used  to  predict  the  current-den¬ 
sity-vs  .  -time  curves  at  the  other  ion  concentrations , 
Figure  15  shews  the  comparison  of  the  predicted  and  ex¬ 
perimental  current-density-vs , -time  curves. 

The  reaction  rate  constants  at  the  applied  overpo 
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Figure  I'Jt  Experimental  and  Computer-Simulated  Curren  t-Densi  ty-vs  .  -T  ime 

Curve:)  with  Various  Nitrate  Ion  Bulk  Concentrations  in  the  Case 
of  Ivl  im:  I  roly  te  Con  taining  No  Cadmium  Ion. 
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tential  of  -0.30V  were  not  pursued.  The  reason, which 
was  discussed  previously,  is  that  at  this  high  applied 
overpotential  more  complicate  sequences  of  reactions 
take  place.  This  simple  reaction  expression  can  not 
describe  the  phenomena. 

The  heterogeneous  rate  expression  so  determined 
is  believed  to  be  a  correct  one.  This  claim  is  supported 
by  the  fact  that  the  rate  constants,  and  Kg,  are 
independent  of  the  bulk  concentration  at  a  given  poten¬ 
tial.  This  is  indeed  the  case  as  was  shewn  in  Figures 
13  and  15. 

6.2  The  Determination  of  the  Homogeneous  Precipitation 
Reaction  Rate  Constant 

Higher  current  density  was  observed  when  cadmium 
nitrate  was  used  instead  of  potassium  nitrate  as  the 
electrolyte  ir.  the  electrochemical  cell.  The  hydroxy 
ion  produced  by  the  reduction  of  nitrate  ion  copreci¬ 
pitated  with  the  cadmium  ion  in  the  solution.  This 
increased  the  reduction  reaction  rate  and  led  to  higher 
current.  This  homogeneous  reaction  was  assumed  to  be 
linearly  dependent  on  the  degree  of  supersaturation  of 
cadmium  hydroxide  (see  Equation  (1^)).  Since  the  beta- 


roger.eous  reaction  rate  constants  X„  and  were  deter¬ 
mined,  the  homogeneous  rate  constant  d  could  be  evalu¬ 
ated  'ey  fitting  the  second  set  of  experimental  data 
with  the  theoretical  model  to  determine  the  value  of 

Figure  16  shows  the  current— density-vs . -time 
curves  with  various  d  values  while  the  other  parameters 
were  held  at  a  constant  value.  The  current-density-vs . - 
time  data  obtained  in  a  0.0025M  cadmium  nitrate  solution 
at  the  applied  overpotential  of  -0.40V  is  also  plotted  in 
Figure  16.  The  reaction  rate  constant,  > ,  was  then  eva¬ 
luated  approximately  from  Figure  lc.  The  most  probable 
homogeneous  reaction  rate  constant  was  determined  to  be 
50  sec-1. 


k=  100 


ity-vt; . -Time  Curves  with  Fixed 
ue  Reaction  Rate  Coikj tault;. 
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CHAPTER  7 

DISCUSSIONS  AND  CONCLUSIONS 

7 • 1  Surface  Concentrations  of  Various  Ions  as  Junctions 
of  Time 


The  knowledge  of  the  surface  ion  concentrations  can  pro 
vide  information  concerning  the  electrochemical  deposition 
process.  Although  the  surface  ion  concentrations  were  not 
available  from  the  experiment,  that  information  could  be 
generated  by  computer  simulations. 

Figure  17  shows  the  surface  concentrations  of  NO-"  and 
OH”  ions  as  functions  of  time  in  the  absence  of  cadmium  ion 
in  the  solution  at  two  potentials,  -0.40V  and  -0.60V  from 
equilibrium  potential  and  at  a  N0^~  bulk  concentration  of 
S.C05M.  In  general, the  surface  concentration  of  NO,"  ions 
decreases  as  time  increases,  while  the  surface  concentration 
of  OH"  ions  increases  as  time  increases.  This  is  due  to  the 
electrochemical  reaction  (23),  in  which  one  nitrate  ion 
reacts  with  two  electrons  and  which  produces  three  hydroxy 
ions . 

For  tne  case  of  -O.m-CV  overpotential,  the  CH”  concen¬ 
tration  increases  sharply  from  the  bulk  value  of  nearly  cere 


KEY 
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to  about  0.5  dimensionless  concentration  units  in  less  than  1 
msec  before  the  hydroxy  ions  generated  are  removed  effectively 
by  diffusion.  Thereafter,  the  OH-  concentration  increases  at 
a  slower  rate.  For  ions,  the  surface  concentration  drops 

to  about  0.72  dimensionless  concentration  units  from  the  bulk 
value  of  i.O  in  less  than  1  msec  and  then  decreases  at  a 
slower  rate  when  the  diffusion  can  effectively  supply  the 
reactant  for  the  electrode  reaction.  The  sharp  changes  of  the 
surface  concentrations  of  and  OH  ions  cause  the  snarp 

decrease  of  the  current  density  in  a  short  time  period  in 
the  beginning  of  the  electrolysis. 


At  a  more  negative  overpotential  of  -C.60V,  the  increase 
in  surface  CH-  concentration  and  the  decrease  of  concen¬ 

tration  are  even  more  significant  and  last  over  a  longer  time 
period.  The  surface  concentration  of  NC^  drops  to  about 
0.5  dimensionless  concentration  units  and  OH  concentration 
increases  to  about  0.57  dimensionless  concentration  units  in 
2.0  msecs.  The  behavior  is  the  consequence  of  the  higher 
overall  reaction  rate  at  -0.60V.  The  nitrate  ions  were  con¬ 
sumed  at  a  higher  rate  which  produced  more  hydroxy  ions, 
longer  induction  time  is  needed  for  diffusion  to  effectively 
supply  the  reactant  from  the  bulk  solution  and  remove  tr.e 
oroduct  from  the  electrode  surface. 


Similar  profiles  of 


sun  ac  i 


cor 


an d  OH 
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at  lower  bulk  NO^~  concentrations  -were  also  obtained.  These 
are  shown  in  Figures  18  and  19* 

Mien  the  electrolyte  contains  cadmium  ions,  the  homo¬ 
geneous  precipitation  reaction  of  cadmium  ions  and  hydroxy 
ions  near  the  electrode  surface  promotes  the  electrochemical 
reaction  on  the  electrode  surface  in  the  cathodic  direction. 
This  causes  the  nitrate  ions  to  be  consumed  at  a  higher  rate 
than  the  case  when  no  cadmium  is  present.  This  is  shown  in 
Figure  20  where  the  homogeneous  reaction  rate  constant  is  50 

4 

sec~l.  It  is  interesting  to  note  that  the  cadmium  ion  surfac 
concentration  increases  initially  and  then  starts  to  decrease 
This  phenomenon  is  explained  as  follows.  Immediately  after 
the  electrolysis  began,  N0^~  ions  were  consumed  by  the  elec¬ 
trochemical  reaction  which  produced  an  N0^“  concentration 

gradient  made  the  NO^"  ions  in  the  bulk  solution  diffused 

+2 

toward  the  electrode  surface.  The  Gd  ions,  on  the  other 

hand,  moved  with  the  N0^~  ions  in  the  same  direction  in  order 

to  maintain  the  elecroneutrality  condition.  In  the  meantime, 

the  OH  concentration  was,  however,  not  so  high  as  to  consume 
+2 

all  the  Cd  ions  brought  in  by  the  transport  process.  The 
cadmium  ions  was  thus  accumulated  and  resulted  in  a  higher 
value  than  its  bulk  condition.  As  scon  as  the  surface  concern 
tration  of  CH  ion  raised  to  some  value,  the  surface  concen¬ 
tration  of  the  Cd  ^  ions  starts  to  decrease  because  of  the 
abundant  supply  of  OH"  ions  consuming  the  cadmium  ions  at 
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this  time. 


For  the  nitrate  ions,  the  surface  concentration  also 
decreases  as  time  increases  at  a  rate  somewhat  larger  than 
that  for  the  case  of  solution  containing  no  cadmium  ion.  This 
is  shown  in  Figure  21.  Figure  21  also  shows  the  difference  in 
the  surface  concentration  of  the  hydroxy  ion  between  those  two 
cases,  i.e.,  in  the  presence  and  in  the  absence  of  cadmium 
ions.  The  initial  behaviors  of  two  cases  are  very  similar. 
After  this  induction  period,  the  ion  concentration  for  the 
case  of  a  solution  which  contains  cadmium  ions  decreases  at 
a  more  gradual  rate  as  a  result  of  the  homogeneous  precipi¬ 
tation  reaction  between  cadmium  and  hydroxy  ions. 

7-2  Concentration  Profiles  of  Various  Ions  at  Various  Time 


For  the  case  of  solution  containing  no  cadmium  ion,  the 
nitrate  ions  were  consumed  and  the  hydroxy  ions  were  produced 
on  the  electrode,  the  concentration  gradients  were  established 
which  made  the  nitrate  ions  diffused  from  the  bulk  solution 
toward  the  electrode  surface  and  the  hydroxy  ions  diffused 
away  from  the  electrode  surface  toward  the  bulk  solution. 
Figures  22  through  27  show  the  concentration  profiles  of 
and  OH  ions  with  different  applied  overpoter.tials  at  some 
selected  times.  From  these  figures,  one  can  see  that  the 
diffusion  thickness, £  ,  increases  with  time.  The  diffusion 
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Figure  22  Concentration  Profiles  of  Nitrate  and  Hydroxy  Ions  at 
Various  Time  for  the  Case  of  Solution  Containing  No 
Cadmium  Ion. 
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Figure  24s  Concentration  Profiles  of  Nitrate  and  Hydroxy  Ions  a 
Various  Times  for  the  Case  of  Solution  Containing  No 
Cadmium  Ion. 
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Various  Times  for  the  Case  of  Solution  Containing  No 
Cadmium  Ion. 
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centration  Profiles  of  Nitrate  and  Hydroxy  Ions  at 
ious  Times  for  the  Case  of  Solution  Containing 


thicknesses  at  various  time  are  listed  in  Tables  2  and  3> 

The  concentration  profiles  of  N0^“,  0H~  and  Ca i’2  ions 
at  various  time  for  the  case  of  solution  containing  cadmium 
ions  are  shown  in  Figures  28  and  29.  The  shape  of  NO  ~ 
concentration  profiles  is  the  same  as  that  for  the  case  of 
containing  no  cadmium  ion  but  the  N  0  ^  ~  concentration  on  the 
surface  is  somewhat  lower  than  that  for  the  case  of  contain¬ 
ing  no  cadmium  ion.  The  diffusion  thickness  increases  as  time 
increases , 

In  the  previous  discussion  of  the  surface  concentration, 
the  surface  concentration  of  0H~  ions  increases  in  the  beginn¬ 
ing.  During  that  time,  OK"  ions  diffuse  from  the  electrode 
surface  to  the  bulk  solution.  Figure  28  shows  th at  the 
surface  concentration  increases  and  the  diffusion  layer  thick¬ 
ness  also  increases  as  time  increases.  Then  the  surface 
concentration  starts  to  drop  and  the  diffusion  layer  thickness 
stays  at  about  the  same  value.  This  is  show  in  Figure  30. 
Table  4  snows  the  diffusion  layer  thickness  for  the  diffusion 
of  OH"  ions  as  a  function  of  time.  In  Table  4,  one  can 
see  that  the  diffusion  layer  thickness  stops  increasing  as  a 
result  of  the  decreasing  surface  concentration  of  Oh’"  ion. 

+2 

The  Cd  concentration  profiles  are  shown  in  Figure  29 > 
As  the  profiles  indicate,  for  each  specific  time,  when  the 
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Table  2:  Diffusion  Thickness  as  a  Function  of  Time 
for  the  Case  of  Solution  Containing  No 
Cadmium  Ion.  The  Applied  Cverpoten'ial  is 
-0.40V.  The  Nitrate-Ion  Bulk  Concentration 
is  Q.005M. 


Time  (sec) 


0.1235  X  10 


-5 


0.6121  X 10 


-5 


0.1758  X10 


-4 


0.3572  *10 


-4 


0.1033  X 10 


-3 


0.4737  X10~3 


0.1258  XlO 


-2 


0.7417  X10' 


0.3946  XlO 


-1 


0.7032  *10 


-1 


Diffusion  Thickness( cm) 


0.25  XlO 


-4 


0.45  X10"4 


0.80  X 10 


—*4' 


0.15  x 10 


0.25  x 10 


-3 


0.50  x 10"3 


0.80  XlO' 


0.20  XlO 


1  o-2 


0.45  X 10 


-2 


0.60  X 10 


-2 


90 


Table  3:  Diffusion  Thickness  as  a  Function  of  Time 
for  the  Case  of  Solution  Containing  No 
Cadmium  Ion.  The  Applied  Overpotential  Is 
-0.60V.  The  Nitrate-Ion  Bulk  Concentration 
Is  0.005M. 


Time  (sec) 


0.1313  XlO"5 


0.4982  X 10 


-b 


0.1642  X10 


zir 


0.3318  XlO' 


0.1298  Xl0~3 


0.4407  X10~3 


0.1258  X10 


-2 


0.8718  X10 


-2 


0.3229  X10 


-1 


0.7004  X 10' 


Diffusion  Thickness  (cm) 


0.25X10 


-4 


0,40  X 10 


-4 


0.70  X 10' 


0.10  x 10' 


0.25  X10 


0.45  XIO"3 


0.75  X10"3 


0.25  Xio 


-2 


0.45  XIO' 


0.65  Xio 


-2 
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30 i  Concentration  Profiles  of  Hydroxy  Ion  When  the  Surface 
Concentration  of  Hydroxy  Ion  Starts  to  Drop. 
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Table  4;  Diffusion  Thickness  of  0K~  Ion  as  a  Function 
of  Time  for  the  Case  of  Solution  Containing 
Cadmium  Ions.  The  Applied  Overpotential  Is 
-0.40V.  The  Nitrate-Ion  Bulk  Concentration  Is 
0.005M.  The  Cadmium-Ion  Bulk  Concentration  Is 
0.0025M. 


Time  (sec) 

Diffusion  Thickness ( cm) 

0.1236  X10“3 

0.25  X10"4 

0.6128  X10-5 

0.45  X10'4 

0.1798  X10-4 

0.80  X10“4 

0.3513  XlO'4 

0.15  X10"3 

0.1017  X io~3 

0.25  X10-3 

0.4856  X 10"3 

0.50  X 10'3 

0.1148  X10-2 

0.70  X10-3 

1  0.8811  XIO"2 

0.85  X10_J 

0.4 050  X 10"1 

0.85  xio'3 

0.7031  X10-1 

0.35  X 10“3 
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concentration  near  the  electrode  is  greater  than  the  bulk 

+2 

concentration  of  Cd  ,  its  value  versus  distance  from  the 
electrode  decreases  slowly  in  the  beginning,  then  faster,  and 
more  slowly  again  before  finally  reaching  its  bulk  value.  At 
the  time  when  the  concentration  values  near  the  electrode 
surface  are  still  larger  than  the  bulk  concentration,  the 
concentration  decreases  to  some  minimum  value  at  some  distance 
from  the  electrode  surface  and  then  increases  until  it  reaches 
the  bulk  value.  The  location  of  the  minimum  concentration 
point  moves  toward  the  electrode  surface  as  time  increases. 

At  sufficiently  long  time,  the  minimum  point  finally  reaches 
the  electrode  surface. 

7 . 3  Conclusions 


The  transport  model  for  the  electrode  kinetics  and  the 
homogeneous  reaction  presented  previously  has  already  explained 
very  successfully  the  characterestics  of  the  current  density  as 
a  function  of  time  for  both  the  case  of  solution  containing 
cadmium  ion  and  that  without  cadmium  ion.  This  model  was  then 
used  to  simulate  the  concentration  profiles  for  various  ions 
involved  in  the  reactions.  The  behavior  of  the  various  ions  on 
the  electrode  surface  and  in  the  solution  was  obtained.  The  ii 
fusion  processes  were  also  known  after  examining  these  concen¬ 
tration  profiles  carefully.  These  can  be  summarised  as  follows 


1.  The  current  density  decays  with  time  after  a  double-layer 
charging.  The  time  for  this  charging  is  so  short  that  we 
can  not  observe  it  by  the  experiment.  The  decay  for  she 
current  density  is  a  result  of  the  depletion  of  the  nitrate 
ions  on  the  electrode  surface. 

2.  The  current  density  for  the  case  of  a  solution  containing 
cadmium  ions  is  greater  than  that  for  the  case  of  solution 
containing  no  cadmium  ion.  This  is  because  the  precipitaticr 
of  Cd(CH)2  promotes  the  electrochemical  reaction  or.  the  elect 
rode  surface  by  consuming  OK”  generated  in  this  reaction. 

3.  The  current  density  depends  on  the  bulk  concentration  of 


nitrate  ions  and  the  applied  overpotential.  Increasing 
the  NO^-  bulk  concentration  or  the  applied  overpotential 
will  increase  the  current  density. 

4.  The  applied  overpotential  is  one  of  the  factors  which 
changes  the  reaction  rate  of  the  electrochemical  reaction 
on  the  electrode  surface.  Increasing  the  applied  over¬ 
potential  will  increase  the  overall  reaction  rate. 

5.  For  the  case  of  solution  containing  no  cadmium  ion,  the 
surface  concentration  of  N0^~  decreases  as  time  increases 
while  the  surface  concentration  of  OH"  increases  as  mime 
increases.  This  is  because  NC^“  ions  are  consumed  to  produce 
the  OH-  ions  during  the  electrochemical  reaction  on  the 
electrode  surface. 

6.  For  the  case  of  solution  containing  cadmium  ion,  the 

surface  concentration  of  icn  also  decreases.  The 

decreasing  rate  is,  however,  faster  than  chat  for  the  c--- 
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oi'  solution  containing  no  cadmium  ion.  The  surface 
concentration  of  OH  ion  increases  in  the  beginning  arc 
reaches  a  maximum.  When  the  production  rate  of  CH~  due  to 
the  electrochemical  reaction  is  less  than  the  consumption 
rate  of  0H“  due  to  the  precipitation  reaction,  the  surface 
concentration  of  CK_  starts  to  decrease.  To  maintain  elec- 
troneutrality ,  the  surface  concentration  of  Cdr^  ions  v/iil 
exceed  its  bulk  value  and  then  decrease  to  the  value  lower 
than  its  bulk  value. 

7.  The  diffusion  processes  for  the  case  of  solution  containing 
no  cadmium  ion  is  that  the  nitrate  ions  diffuse  from  the 
bulk  solution  toward  the  electrode  surface  while  the 
hydroxy  ions  diffuse  in  the  opposite  direction. 

8.  For  the  case  of  solution  containing  cadmium  ions,  the 
diffusion  directions  for  both  the  nitrate  ions  and  hydroxy 
ions  are  the  same  as  that  for  the  case  of  solution  containing 
no  cadmium  ion.  However,  the  cadmium  ions  diffuse  toward 
the  bulk  solution  in  the  beginning  and  then  finally  change 
direction  toward  the  electrode  surface. 

7  •  4  Recommendation  far  Future  Works 

Although  the  model  has  successfully  explained  the  cases 
of  -0.407  and  -O.cQV,  it  is  not  likely  valid  for  the  case  cf 
-0 .  oC 7  due  to  the  unexpected  characcerescics  of  the  currer.t- 
de.nsity-vs . -time  curves  and  the  comparison  of  the  current 


density  with  the  limiting  case.  V/hen  the  applied  cverpoter.tial 
is  as  high  as  -0.80V,  the  reaction  mechanisms  and  electrode 
behavior  should  be  identified  by  doing  some  further  experiments 
Furthermore,  the  current-density-vs . -time  data  at  -C.cOV  applie 
overpotential  for  the  case  of  solution  containing  cadmium  ions 
is  not  available  yet.  It  can  be  obtained  by  repeating  the 
experimental  work  we  decribed  before.  Once  the  data  is  obtained 
the  model  can  be  used  for  comparison  with  the  experimental 
data. 


Finally,  the  configuration  of  the  electrode  in  the 
impregnation  process  is  not  as  simple  as  we  have  used  in  derivi 
the  model.  It  is,  in  fact,  a  porous  electrode  in  the  real 
case.  If  a  micro-porous  electrode  is  possible  to  obtain, 
the  experiment  can  be  repeated  to  obtain  the  current-density- 
vs. -time  data.  Cnee  the  data  are  available,  the  transport 
processes  can  be  examined  for  the  condition  of  porous  electrode 
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APPENDIX  A.  1 


THE  "BASIC"  MAIN  PROGRAM  USED  IN  EXPERIMENTAL  WORKS 
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oo 

35 

63 


13  F'EM-»- 
13  REM* 

_0  i-'EM*  iHIS'  I S  THE  £:RS I C  PROGRAM  USED  IN  THE  C LFPEHT  '6. 

63  P.EM*  TIME  TRRNSIE'NT  E'.-'.FERIMEHT .  THE  "MUG"  IN  ET.-iTEMENT 
30  REM*  L6u  IS  A  SAMPLING  SUBROUTINE  URITTEN  IN  ASSEMBLY 
33  REM*  LANGUAGE .  WHEN  "MUG"  IS  EXECUTED-  THE  CPU  STARTS  TO 
40  F'EM*  SAMPLE  THE  CURRENT  OUTPUT  AT  EUEF.Y  CONSTANT  T I  ME 
43  REM*  INTERNAL.  THE  ANALOG  SIGNAL  IS  IMMEDIATELY  CONCERTED 
REM:*  INTO  A  DIGITAL  DATA  BY  A  ID  EOT  ADC  AND  THE?''  STORED 
r'EM*  IN  i HE  MEMORY .  UHEN  THE  NUMBER'  UF  TIMES  IT  SAMF’LES 
REM*  PEACHES  A  PRESET  UALUE .•  THE  SAMPLING  ROUTINE  IS 


6-6 

REM:* 

DEC:  IMAL 

NUMBER'S  AND  RR.  INTZD  OUT  E;V  A  LI  NE 

PR I N i HR 

C1 

REM;* 

THE  UR 

R  IPlBLES  USED  I?l  TH I S  PROGRAM  ARE 

7*0 

REM-*: 

72 

REM* 

o> 

— LOOK  SE  •  i  I  No 

r  4 

REM* 

'  ' ,  1  , 

=SnMFL I NU  TIME  I iN  i  - ;-1 -'hL  .N  MIi_.Rus 

ECONO 

“*  _r 

REM:* 

>  ‘ .  2  •' 

=UHANNEL  NUMBER  INHERE  THE  U UP  RENT 

x  ED 

r  o 

REM* 

'SAMPLED  PROM 

so 

REM  * 

M 

=NUMZ E.R  OF  POINTS  TO  GE  SAMPLED 

•  -■2 

REM:* 

THE  AP 

GUMENTS  IN  "MUG"  SUE-ROUTINE  ARE 

94 

REM:* 

96 

REM;* 

IP  S 

=0  THEN  EXECUTE  TEST  PROGRAM 

•Z’Zf 

REM:* 

IP  s 

=  1  THEM  GO  TO  MUGGER 

90 

REM* 

N=NI. 

MBER  OP  UALUZS  PASSED  Iri  APR  AY  ?v. 

0  1 

92 

REM* 

M=NL 

ME'EP:  OP  ! vALUES  TO  EE  RETURNED  EM 

r  7  7nY  V'. 

i 

REM:* 

IS  THE  INPUT  ARRAY 

94 

REM:* 

?n::o> 

IS  THE  INF'UT  AP  RAY 

96 

REM* 

v  <  >3  y 

IS  THE  OUTPUT  CURRENT  ARRAY 

97 

REM:* 

.30 

REM* 

9’? 

REM* 

t . .T  .... 

1  TO  CGf  <F  I G 1  SCREE?  ?  1 

1 3u  D I M  y'.  <  1 0  .•  V 3000  >  >  Z 3000'  '*  .•  C 300  1  .■  T 1  30* 
1  ?0  S— 0 
ZOO  N=3 
GO  1  P  EM'* 


GOG 

REM* 

204 

REM*:**-*:*  INF'UT 

Z-:  nTn  r-.J-Nr.-r-v 

:G03 

REM* 

GOO 

REM* 

2  i  0 

PRINT  "CLOCK  SET . 

1  .•  " -'f-irir L  I No  r.?L  .NT  w,-- -  PL' 

ZGO 

-  — 

INF'UT  0 • ; ■, \  1  > 

~  _• 

2-0 

~  -  IN  "uHANNEL  iNUi 

■iBER .■  "  •  "NUMBER  CP  i C I :  «TS  •  "  ‘S 

G30 

INF'UT  M1  G.>  ■  i?..  P 

G31 

REM* 

02 


4—  ! 

234 

CFLL  7HE  MUG  SUBROUTINE  PND  S7nR~  70  SnriPLt.  "*  ** 

233 

REM* 

236 

REM* 

2  60 

ooit.mug. 

S  >  N  .•  '■  0  ’>  .■  M  .•  V  <■.  0  ’>  l1 

261 

REM* 

REM* 

264 

REM****  r 

GNUER7  7HE  CURR.EH7  DFTri  IN7G  DcGIMPL  NUMBERS  <••"*" 

263 

PEM* 

REM* 

»  'J 

FOP  1=073 

M- 1 

•x-CnJ 

Zv  i  a 

.>  —  15  >  •••'  1 6  ••■'2043 

200 

next  i 

REM* 

_  ,"***-4 

PEH-*: 

300 

(>  o'"*** 

CONNECT  7HE  CGMPU7ER  70  THE  LINE  PRINTER  "■-••"* 

302 

Di' -IfliL-1  J  Hu  i  »■'  i  Hu,  2-ijrT  iii.  c.  r  >r  z_n  rtf..*. 

304 

■vfc.M***" 

-J  .013  UM2  1  O'  GEOGnE  "HE  CURRENT  lEriSIi",-' 

303 

REM"*** 

r  r  iNi  OijPr  tH  i  Nttl'lSl  i  V  O'-.  71  ME  DPiTni  *•*•<•" 

306 

REM* 

307 

REM 

310 

CONFIG 77 

320 

PPIH7  7nE 

13.-’  ••  "  '  I ME"  •  7 ,-iB 43 >  .•  "UURRENY  DENS I TV" 

330 

PRIN7  7 Fit' 

v  16.>  .•  1  EEC-  >  "  .•  TPS »l SO .> ..  "  <  FiMP.-1 'CM2  >  " 

340 

FRIN7 i FR I 

NT 

330 

FOR  1=1  - 

J,  i  o' 

3c0 

7  ■'  I  >  =P  10 

I 

370 

C-:  I  '=-2v  I 

-•  1  >  i-r  0.  1  •  0.  0 1 3 

360 

FRIN7  7FiE: 

•.  1 3'1  .•  7'  I  >  ..  7  Ti  E-  <  50'  >  .•  '2  <  I  > 

3  ?0 

NEXT  I 

400 

FOR'  1  =  1  7 

3  M  ,  20'-  1 

410 

7'  I  <=►.'  10 

30-1 30  »• lO:+:  I  > 

420 

C  I  "'=-3'  2 

— '  r  *.  1  1  *  *i "  •  2  ■••  •  k  _■ 

430 

FR I N7  7 FE 

■  *-•  •  .«  k  .'-L''  Z-*J  . '  .•  ■-  x 

440 

ne: '7  I 

430 

END 

! 
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APPENDIX  A. 2 

THE  ASSEMBLY  SUBROUTINE  "MUG"  USED  IN  THE  MAIN  PROGRAM 


Address  Contents  Instructions  Comments 
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ov-v 


Contents  Instructions  Comment.: 


0 C9  RCT  Return  to  main  program 


APPENDIX  3 

LISTING  OF  EXPERIMENTAL  DATA 
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:a  cue 


3.3  Curran* -Density- vs .  -Tine  Tata  for  the  Case  of 
Solution  Containing  Ho  Cadmium  Ion.  .titrate- 
Ion  3ulk  Concentration*  O.QCfM,  Initial  pH=  3.0, 
Aoplied  Overpotential*  -0.30V,  Rest  Potential 
=*  -0.35^. 


CCS 
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Table  3. 


^urrent-C  er.  s  i  ~y  -  vs 


Solution  Containing  ': 


-me  -a:a  ;or  one 
o  Cadmium  Ion.  :iitr 
Ion  Bulk  Concentration5  0.C025M,  Initial  ?H=3 • C , 
Applied  0verootential=*0 . 6CV,  Rest  Potential 
■  -0.35^. 


CiJ 
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Tacxe 


3.6  Curr  =r.t-Der.sity-vs  . -Tine  Tata  for  tr.s  Cas 
Solution  Containing  No  Cainiur.  Ion.  Nitra 
Ion  Bulk  Concentrations  0.2 C2;~A,  Initial 
Aoolied  Overootential=  -C.SCV,  Best  Poten 
=*  -0.35^-V. 


f  • 

d 


m  I  ’  Pi  P 


116 


■a-Ie  3.9  Current-Ber.sity-ys  .-Title  Tata  far  tr.e  rase  of 
Solution  Containing  No  Cadoiiun  Ion.  Nitrat 
Ion  Bulk  Concentration3  C.CC125M,  Initial 
Applied  Overootential3  -O.cCV,  Best  Potential 
=  -0.354V. 


CURRENT  Per  15 1 T’ 

•  rtf'lp  CM2 


.  000 1 
.  0002 
.  0003 
.  0004 
.  0003 
.  0006 


.  WOOF 


.  oolJ6 
.  000'? 


•n 

3 
3 
3 
3 

.  00 1 3 
.0016 
.  00 1 7 
.  00 1 6 
.001'? 


164.043 
141.976 
123.  l‘?6 
1 1 1 . 1 37 
102. 333 


Ktt-1  ,•  M. 


.  004 

33 . 3 1 1 2 

.  006 

3o . 4'i36 

.  003: 

25.3164 

.01 

22. 3116 

.012 

2 1 . o336 

.014 

13. 1536 

.016 

13.4043 

.  0 1 3 

17.2776 

.  02 

16.5263 

•  0  7,  •=;  t- 

! 


•'ll  a. 
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Table  3.12  Background  Current  Density  of  Water.  Initial 
pK=  3.0,  Applied  Overpotential3  -0.90V,  Rest 
Potential3  -0.35^. 


T I  ME  CURRENT  DENS I TV 

CSEO  < m  HMP--  CM2  > 


.  QOO  l 
.  0OO2 
.  0003 
.  0004 
.  0005 
.  0006 
.  0007 
.  OOOS 
.  0009 
.  00 1 
.0011 
.0012 
.  00 1 3 
.  00 1 4 
.0015 
.  00 1 6 
.  00 1 7 
.  00 1 3 
.  00 1 9 
.  002 
.  004 
.  006 
.  003 
.01* 
.012 
.014 
.016 
.013 
.  02 
.  022 
.  024 
.  026 
.  023 
.  03 

.  Uw2 
.  034 
.  036 
.  033 

.04 


247. 145 

135.216 

93.5246 

72. 1153 

60.0961 

52 . 534 

47. 7013 

43.9453 

41.316 

33.6363 

37 . 1 344 

3<=;  ^,32 

34 . 1 796 

33.4234 

32.3016 

31.5504 

30 . 7992 

30 . 043 

29. 6724 

23.9212 

22.536 

19. 9063 

13.0233 

16.5263 

15.3996 

14.6433 

13.3972 

13.5216 

12.7703 

12. 3947 

11.  q.436 

1 1 . 6436 

10.3923 

1 0 . 3923 

10.5167 

10.5167 

9. 76561 

9.  76.561 

ci_  -yq 


cl  O 
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Table  3.1c  Current-Density-v; 


■  line 


:  o: 


L3e  or 


Solution  Containing  Cadmium  Ions.  Nirrare- 
lon  Bulk  Concencration=  0 .0023M( _ Cadmi 
Bulk  Concentration3  0.00125M,  Initial 

Apolied  Overpotential =  -Q.4CV,  Has*  ?o - 

=  -0.35^V. 


■JS  P.P 


O 


126 


Taels  3.19  Current -Density- vs . -T 

Solution  Containing  Cadmium 


ata  tor  mr.e  Case  ci 
ms.  Nitrate- 
Ion  Sulk  Concentration3  0  .CGlSo^i  Cadmium- 1 cr. 
Sulk  Concentration3  0,0006252/1,  Initial  ?H=  3.Q, 
Apolied  Overcotential=  -0.40V,  Rest  potential 
=  -0.354V. 


t-i 


Bulk  Concentration3  0. 00625'*!,  Initial  prl= 
Aopiied  0verootential=  -0.807,  Rest  Potential 
=' -0.354V, 


M 
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The  method  of  weighted  residuals,  KWR,  is  a  general 
method  for  obtaining  solutions  to  equations  of  change,  in 
our  case,  Fick's  second  law.  In  the  MWR ,  one  assumes  a  trial 
function,  usually  a  set  of  weighted  polynomials , substitutes 
this  trial  function  into  the  differential  equation  and  then 
selects  the  coefficients  of  the  polynomial  terms  by  specify¬ 
ing  that  the  residual  be  zero,  on  the  average,  at  certain 
points.  If  one  evaluates  the  differential  equation  at  the 
zeros  of  an  orthogonal  polynomial ,  the  residual  will  of 
necessity  be  exactly  zero  at  these  collocation  points.  3y 
increasing  the  number  of  collocation  points,  the  trial 
function  would  satisfy  the  differential  equation  at  more  and 
more  points^ 1 ^ . 


(1)  3. A.  Finlayson  i  Richard  Bellman  (Id,),  'mathematics 

in  Science  and  Engineering’,  Vol.  8?,  Academic  Press, 
New  York,  1972,  pp  9 . 
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with  the  boundary  conditions  of  the  form 

y(o,t)=g(t)  (3) 

y ( 1 , t  )=h( t )  (4) 

and  initial  condition  such  as 

y(x,0 )=p(x)  (5) 

For  a  mass  transfer  problem  D  would  be  a  dimension¬ 
less  diffusion  coefficient  and  f(y)  would  be  homogeneous 
reaction  term.  Geometric  considerations  indicate  that 

o£  =  £3=0  and  mathematical  considerations  indicate  that  a 

( 2 ) 

suitable  but  not  unique  trial  function  isv  ' 

N 

y(x,t)a(l-x)y(0,t)+xy(l,t)+x(l-x)ZIai(t)P.  .(x)  (6) 

d"  1 

The  ? 4  1  of  equation  (6)  refers  to  one  member  of  a 
0  **  x 

complete  set  of  polynomials.  Table  0.1  lists  some  Legendr 
polynomials ,  the  special  case  of  equation  (1)  -where  ot  =  3 
ar.d  their  roots  .  Note  that  a  Legendre  polynomial  of  degr 
has  N  real,  distinct  roots  or.  the  interval  0<x<  1. 


\Z)  z.k.  Finlay  son  -  Richard  Bellman  'Ed 
in  Science  and  Engineering' ,  V ol,  37, 
New  York,  1971,  pp  105 ■ 


'Mathematics 
idenic  Press, 
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TA3LZ  C.  1 


Legendre  Polynomial  and  Their 


r*.oots 


(3) 


fl 

PM  - 1 

Roots  of  P.. 

iN  ■ 

1 

-l+2x 

0.500000000 

2 

l-6x+6x2 

0.211324865 

0.788675135 

3 

-l+12x-30x2+20x3 

0.112701665 

0.500000000 

0.887298335 

4 

1-20x+90x2-140x3+70x4 

0.069431844 

0.330009478 

0.669990522 

0.930568156 

(3)  John  Villadsen,  'Selected  Approximation  Method  for 
Chemical  Engineering  Problems',  Reproset,  Copenhagen, 
Denmark,  1970,  pp  A3  and  A4. 
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I x  is  evident  that  regardless  of  the  value  of  N  cho¬ 
sen  one  can  always  rewrite  eqn.(6)  in  the  general  form 

M  +  2 

y(x,t)»y  xJ"1b,(t)  (?) 

*-  '  ~  O 

<j  =  i 

One  could  at  this  point,  as  is  normally  the  case  in 

using  MWR,  obtain  the  coefficients  b-(t)  by  substitution  of 

j 

eqn.(?)  in  the  differential  equation  whose  solution  is  de¬ 
sired.  However,  the  computer  programming  is  greatly  simpli¬ 
fied  when  they  are  written  in  terms  of  the  solution  of  the 
differential  equation  at  the  collocation  points,  y(x. ,t) 
where  i=i , 2 , 3 • • .N+2 ,  and  each  x,  is  one  of  the  N  roots  of 
the  polynomial  or  one  of  the  two  boundary  values.  So  in 
terms  of  the  N+2  collocation  points  we  may  write: 

N+2  . 

y(x, ,t)=  V xP_1b  .(t)  (8) 

U 

We  now  have  an  approximating  polynomial  function  which 
we  will  force  to  fit  our  partial  differential  eqn.(2)  at 
the  N+2  collocation  points  (i.e,  cnoose  the  b,(t)  so  that 

J 

our  partial  differential  equation  is  satisfied).  In  order 
to  do  this  we  will  need  both  the  first  and  second  derivatives 
witn  respect  to  x.  These  derivatives  can  be  calculated  ex¬ 
plicitly  in  terms  of  x  since  we  have  already  defined  the 
dependence  of  y  upon  the  spatial  coordinate ,  ecr. .(  o ;  and  (8). 

It  is  easily  seen  that  if  one  considers  ail  the  co- 


1 3^ 


liocation  points,  x^,  equation  (5)  represents  one  xen'cer  o 
a  set  of  equations: 


N+2 

y(x1 , t)  =  y~ xj~1o j( t)  (9) 

]=' 


Cne  can  see  that  the  left-hand  side  is  a  vector  of 

the  trial  solution  evaluated  at  the  collocation  points. 

•  « 

Since  the  term  involving  depends  on  both  i  and  a 

matrix  will  represent  this  term.  The  b;(t)  will  naturally 

O 

be  a  vector.  If  v/e  let  a  single  bar  represent  a  vector  and 
a  double  bar  a  matrix,  then  the  set  of  ears. (9)  through{li 
can  be  'written  as  one  equation  in  this  simplified  notation 

y  ( t / =Q • b( t )  lily 

wnere  Q; 

^  •  o 

The  first  and  second  derivatives  with  respect  to  x 
may  then  be  obtained  from  equation  ill): 
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—  -bR 

dx“ 


(1*0 


where  C 


dxJ_1  i  ,  n  d2xJ_1 

,d  d  i*i  10  ax4  jx. 


it  is  now  simple  to  soive  lor  o^tj  tne  vector  ox  o ..  (t; 

xJ 

coefficients,  since  simple  rules  of  matrix  algebra  apply  in 
this  notation.  Equation  (12)  yields 


o  ( t )  =  Q  "  y  ( t ) 


Substituting  this  into  eqns.(13)  and (14)  gives 


dy(_ 


dx 


Q'x  y(t)=  A  y(t) 


( 15) 


b  y  (  V  )  -  -V  .-1  _  1  ,,  (  X.  \  _  3  .  -  (  —  ; 

dx4 


(16) 


where  A=  C  Q-1  and  3=  D  Q-^  or  in  terms  of  the  i-th  collo¬ 
cation  point,  after  some  simple  matrix  multiplication, 


dy ( x ,  t )  _ 


rt*2 


dx 


=  y  :  i  y  V  x  . ,  t ) 


1 J 


,2.  - . 

—  3 ^  1  j -1.  =  T~ 3  --,-c 

.4  >  “i,J 


dx 


J=l 
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One  can  therefore  express  the  spatial  derivatives  in 
terms  of  the  trial  function  at  the  .'i+2  collocation  points  . 

Equation  (16)  can  be  directly  applied  tc  some  diffu¬ 
sion  process  such  as  that  described  by  ecns.(2)  through  (5). 
Equation  (2)  becomes 


dy  ( t )  _ 
dt 


f(y( t) ) 


(19) 


where  the  vector  dy(t)/bt  represents  the  set  of  time  deri¬ 
vatives  at  the  Nr2  collocation  points,  be  note  that  there 
is  no  explicit  spatial  (x)  dependence  in  eqn.(l?).  3y  assum¬ 
ing  an  explicit  polynomial  in  x,  eqn.(6),  and  evaluating  it 
at  specific  values,  the  collocation  points,  the  problem  is 
reduced  from  a  partial  differential  equation  with  two  inde¬ 
pendent  variables  to  a  set  of  ordinary  differential  equa¬ 
tions  . 


One  can  expand  eqn.(19)  into  its  separate  members  to 
solve  the  problem.  For  the  i-th  collocation  point  we  can 
•write,  according  to  eqr..(l9), 


for 


1  »  ^  »  3  I  •  •  •h",‘2. 


13? 


t 

i 

i 


i 


One  now  considers  the  conditions  at  the  boundaries 
x^=0  and  x^^l.  Since  these  values  are  known  from  the 
original  problem  i.e.  eqns.(3)  and  (4)  one  can  subtituxe 
the  boundary  conditions  in  eqn.(20)  to  obtain 


d7dt’  ~  =  kLBi,l*s(t)"  3i,N-2*h(t)' 


Pi,jy(xj*t)]"  f(y(xi’t)} 

Ja| 


for  i  =  2,3,4. ,.N+1. 


One  should  note  that  g(t)  and  h(t)  are  included  so 
that  the  boundary  conditions  may  be  time  dependent  in  which 
case  g(t)  and  h(t)  would  be  known  or  easily  obtained.  The 
diffusion  problem  has  been  transformed  into  a  sex  of  N 
simultaneous  first  order  ordinary  differential  equations 
in  N  unknowns  with  initial  conditions  and  can  be  solved 
numerically  by  a  number  of  standard  methods  on  a  computer. 

Evaluation  of  the  current  at  an  electrode,  i=  nFAD|^-!  , 

Sx  jx=0 

can  be  obtained  directly  from  eqn.(17)  once  the  concentra¬ 
tion  profile  is  obtained. 


I 


APPENDIX  D 
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PROGRAMS  USING  ORTHOGONAL  COLLOCATION 


METKC I 


1 
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j-*********«*******»*tt**  *«***************«*?  «*«****4****w*-*«««*** 

¥  'J 

C 

c 

THIS  PROGRAM  SOLVES  EQUATIONS  (50)  THROUGH  (55)  n'li.i  Tao 

X 

G 

INITIAL  0  0  ii  D  I  TIG  N  S  USING  SPLINE  COLLOCATION  METHOD  FO  A  I  :i  r, 

c 

c 

CASE  OF  THE  ELECTROL YT  E  CONTAINING  NO  CADMIUM  ION. 

z 

-  c 

EQUATION S  (DO)  AND  (51)  A  HE  DISCRETIZED  FI AST -03 DEE 

C 

c 

OB  D I N  A  3  7  DIFF  EB  ENT  I A  L  EQUATIONS  «HICH  ABE  SOLVED  BY  AN 

w 

c 

SEMI-IMPLICIT  INTEGRATION  SUBROUTINE  CALLED  DIFSUB  USING 

c 

c 

GEAR'S  ROUTINE.  EQUATIONS  (52)  THROUGH  (53)  ABE  SOLVED 

SIJULI  ANS3U3L  i  BY  A  SUBROUTINE  CALLED  0  SY  3  T  M  Nil  ICH  COMES 

•g 

V- 

FROM  I  NS  4.  S  J  jSOU  TIN  E  PACKAGE  u  OBTAIN  THE  SUE  FACE 

c 

0 

CO  NC  EN  T  B  A I  .  0  si  o  0:  NITRATE  AND  HYDROaY  IONS.  EQUATIONS  (5-tj 

c 

J 

AND  (5  5)  ASi  Tn SIB  BULK  VALUES.  THE  CO N C E N 7 Ii AT 1 0 N  P!iOrI..SS 

J 

c 

THUG  03TA1NLD  ABE  THEN  USED  TO  CALCULATE  THE  CURRENT 

c 

c 

DENSITY  jY  :t  UA  TIG  N  (5b).  TiiE  CURES..  T  DENSITY  DATA  ARE 

c. 

G 

THEN  USED  YO  FIT  KITH  THE  EXP  ERIdENT AL  DATA  TO  IDENTIFY 

G 

G 

THE  CATHODIC  AND  ANODIC  REACTION  BATE  CONSTANTS  0?  THE 

c 

r* 

SURFACE  REACTION. 

c. 

c 

V. 

c 

-  NO HEN CL A. USE  - 

- 

>>- 

n.(=  number  of  points  in  x  di..ectio.i. 

N=NrJd3ER  Or'  INTERIOR  COLLOCATION  POINTS. 

c 

r* 

NT= NUMBER  OF  TOTAL  0  CL LOCATION  POINTS  INCLUDING  CNJ 

c 

c 

JCUtJOfioLl  £ *J x T S • 

n 

v- 

Y  ( 1  ,  1 )  ,  Y  ( 1  ,  L)  ,  .  .  .  Y  ( 1  ,  N)  =  CON  CENTS  A  TIO  NS  OF  NITRATE  ION  AT 

G 

r 

S  I  N  TESIOK  COLLOCATION  POINTS. 

c 

v» 

Y  ( 1 ,  N»  1 )  ,Y  (  1 ,  N  2 )  , .  .  .Y  (1 , 2*N)  =CONCEN  TR  AT  10  NS  OF  HYDRO  At 

c 

c 

ION  AT  N  INTERIOR  COLLOCATION 

c 

u 

POINTS. 

V- 

c 

Cl  (1)  ,  C  1  (2)  ,..o1  (NX)  CONCENTRATION  OF  NITRATE  ION  AT  EACH 

c 

o 

POINT  IN  X  DIRECTION. 

c 

w 

C2  (  1)  ,~2  (I)  ,  .  .02  (NX)  CONCENTRATION  OF  H  YDROXY  ION  AT  tAC.i 

c 

G 

POINT  IN  I  DIRECTION. 

c 

\  z 

X S  ( 1 )  3  3U  A  F.iC  t.  CO NC £ N TRATI CN  OF  NITRATE  tON. 

c 

c 

X3  (0  )  =  SU  R  F  AC  -  CON’C  ENTP  AT  ION  or  HYJaOXY  ION. 

G 

!  c 

C 1 3=  BULK  CONCENTRATION  OF  NITRATE  ION. 

G 

t 

C2  3=  3li  LX  C  J  N  C  EN  TRATI  CN  OF  HYD30CY  ION. 

c 

D'1AX  =  D  IF  FUSION  THICKNESS,  CM. 

C 

1 

03  =  LOC  AI 10  N  or  THE  SPLINE  POINT,  G  <  05  <  1  . 

C 

1  c 

DO  3  =  INCREMENT  OF  THE  SPLINE  POINT,  OS. 

c 

II MAX= MAXIMUM  .IME  TC  BE  IN  I  OCR  AT  ED ,  SEC. 

c 

JT  IME=  T I. Ir,  INOaEMENT,  SEC. 

w 

G 

TAJMAX=OIMO.I5  IoNLESS  MAXIMUM  TIME  To  BE  INTEGRATED. 

,'3 

J 

A  rt  i*3i,JU*lSxJlt  MU  O  O  A  L  '1  W 

c 

-  .  J  —  J  1 .1 M  ,i  J  2  sj  »«  L  Z Zi  S  I  1  i\  Z  I  !IC  ^  Z  A  G  j  I 

G 

z 

oA-Ia-a  iu.i  3  ")  Z  3  C  i'  ."i  L  D  2  J .(  i  1 0  G  m 

G 

1 

U  -T  A=r  iACT  :j.i  ur  J"?  OF  NITRATE  ION. 

f 

Jo  D  wdot  IQJ  t\  A  1  u  wJ  .1  O  Ta  .i  k  X.*  M.i  J  Jis, 

G 

Z  ±  ti  ws,  »  IC  ... 

- 

CO.i  1  =  H  oT  oN  D  G  o  -i  -0  US  REACTION  RATE  CJNJTA.iT  IN  CATHJU.C 

"■* 

Z  A  *v  wCa  iw  .«  « 

0  D-CUS  RENT  DE.iSTTY,  AM  P/C  M*  «  z  . 

-i 

- .  • 
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C  THE  3IFFJS1JN  COEFFICIENTS  ■>:  Ini  NITRATE  A.<D  HI ORCaE 

IONS  ARE  1. <#022-5  A  N  0  5.  2  60  E-5C3  I/S  EC ,  3  ESP  ZCTI V  E  „  { .  C 

C  C 

q  turn******  ********  **  ******************************************  *  C 


IMPLICIT  REAL'S  (A-H,  0-2) 

DIMENSION  A  (20,10)  ,0  {20  ,20)  ,C  1  (5  0)  ,C2  (60)  ,  C  (2  0)  ,  I  (8 ,20)  ,  1 2 AX  (2  0)  , 

1  OIF  1  (20)  ,01  E2  (2  0)  ,0  IF3  (2  0)  , ROOT  (20)  , 7ECT  (20)  ,  ERROR  (20)  , 

2  S  A/ 2 (24  ,  20)  ,  ?*  (40  0)  ,  XS  (2)  ,  AX(  1)  ,  WA  (  J)  ,  3AS  (1) 

EXTERNAL  A  J  Xe 

CC.liIGN/LI/CJIE  1 , C j  ZF5 

COi'Ei  ON/L2/OOEF  2 

C  Gild  UN  /  L3  /  N  ,  N  r 

CLKLIGN/L-i/A  ,  5 

COilNON/L5/\S 

CONXOtN  /Lo/COiF  4 

C  U  aM  G  N  /  L  7  /  N  S I  3 ,  N  2  Z  ,  Z  9 

C 0.1(1  ON  /  Ld/C  02”  3, CONI  ,C0  N2 , 3  AO  A  ,  JS I A 

cg;i::cn/lo/y 

-  INPJ,  DATA  - 


c 


READ  (5,3  10) 
R EnJ (5,320) 
READ  (5,33  0) 
HEAD  (5,340) 


l'i  0  ,  :i,  NO,  N  1  ,  AL  ,  3  E,  EPS,  TIE  A  a, 1C, DOS 
C 1  3  ,  C  2 3 ,  0(1  A  X ,  0  T  a.1  2  ,  N  v ,  2 S 
vjAi*i*»,  B  IT  A,  C  C  il  2  ,  CO  N  1 

viSIO,NtE,^P 


- CALCULATE  SG.1E  CONSTANTS  TO  dZ  USED  IN  TNE 


PR  OS  RAO - 


COEF1  =  (  1,  0  02  0—  5 )  /  5-  2  b  0  -  5 
COi ; 2  =  C^  3/  C  1  5 

COE  ?3= ( 1 .  J  020- 5)  *C1 S/D  NA X/ES 
COE F4  =  2  .  *(1.9023-5)  *C1  3*Sb48  7„/0(1A  a/ZS 
C02F5=1./<,3/23 
Da=Di'1Aa/  (Na-1) 

DT  A0  =  (  D.  2  o  J  -  5)  *DTI  XZ/0.1  A  a/0  .1A  a 

T  A 0 0  A  X  =  ( 5 .  2  b  u~5)  *T  Id  AX  /  0*1  A  X/  0.1  AX 

SDTA J= jTA  J 

DETA= 1 .  /  (Na-1) 

nt=n  ♦no*:;  i 

NE2=2*N 


SET  INITIAL  DI 0  ENS CONFESS  TIN Z  =  0. - 


TAJ3 0. 


A  .i  J  .UNK1J  .1  01(1  ENS  I J  N  LESu 
a  j  d  E  0  IN  T . ;  £  IL-js.i.  I  ^  N 

OTAO  .11  =  .  0  o  1  *  L  T  A  U 
DIrtJ.1A=2.0*DTAJ 
«  S I T  2  ( u  ,  j  j  0  ) 


■-  D  E  E I «  E  7  ii  L  :1AaI;VJ.1 
,  i  .1  w  INC  ->  .*£  .»  l  l  3  c 
SUEa^UiiNa  D.roOS  — 


1 


n  n  r,  r,  n  n  n  c,  n  n  r,  n  Ci  n  n  r.  n  n  n 
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- E  V  A  LJ  AT  S  THE  COLLOCATION  POINTS  (BOOTS )  OF  THE 

JACOol  POLYNOMIAL  C?  ORDER  N,  AS  N  ELL  AS  Tel 
FI  a ST  AND  SECOND  DE3I V  AT  IV 25  OF  THE  POLYNOMIAL 
h T  TaE  HOOTS.  - 


CALL  J C 0 0 1 ( N  D  ,  M  ,  N  0 ,  N 1 , A  L , 3  S ,  D I F 1 , 0 1 F2 , DIF3, HOOT) 

WHITE  (6 ,350)  (3001(01,0=1,31) 

- CALCULATE  THE  DISCRETIZATION  COEFFICIENT  MATH  1 1  A - 

1 0  =  1 

DO  20  1=1, NT 

CALL  DFO?  d  (ND,  N  ,  NO  ,N  1 ,1 ,  ID,  DIF  1 ,  DIF  2,  DIF  3,  HOOT  ,  VECT) 

DO  10  0=  1 , NT 
10  A  (I,  J)  =7ECT  (J) 

20  CONTINUE 

-  CALCULATE  THE  DISCRETIZATION  COEFFICIENT  caZaIk  3  - 

ID=2 

DO  4  0  I  =  1  ,  NT 

CALL  D  F  3  ?  a ( N  D , 3 , NO , N 1  ,I,I0,DI?1 ,  DIF2  ,  DIF3  ,  ROOT  ,  V  EC.  ) 

DO  30  0=1, NT 
30  3 (I, 0) = VECT (0) 

40  CONTINUE 

-  XNPJT  THE  VALUES  OF  ARGUMENTS  TO  5E  USED  IN  THE 

INTEGRATION  SUBROUTINE  DIFSU3  - 

Mr  =  1 

OSTA3I=0 
M AID  SR  =  7 
DO  4  5  1=1, N  SB 
7  M  A  I ( I )  =1.0 
45  CONTINUE 
1  =  0 

-  INPUT  THE  INITIAL  CONDITIONS - 

DO  5  0  0= 1 ,  N 
I  0  ,0)  =  1. 

Y  ( 1  , 0  ♦N)  =CO£F2 
50  CONTINUE 
XS(  1)  =  1. 
a  j  U  )  —  vO  .  F  ** 

W H II E ( o , 3*0) 

SO  1 1 .12  =  1 A  U  *  D .»  A*  *  D  U  A  X  /  (5.  2  60-5) 

AND  PRINT  THE  CU  R  3  E  N  T  DENSITY  AT  Tills  .IMS 
CONCENTRATIONS  OF  NITRATE  AND  sYZROIY  lo  ..5 
COLLOCATION  POINTS - - - ■ 

CALL  C  JRENT  ^CD,  A,  Y  ) 


™  Cti.cjLA.il 

Pal N.  THE 
a  l  .HE  N  T 
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L 

c 

c 


*RxTE (o,3o0)  T I  U  E  ,  C  0 ,  X  J  (  1 )  ,  X  a  ( 2 }  ,  X/L  A  0,1123  N,f,u 
*  SITE  (5,700)  (7  (1,1)  ,  L=  1  ,  .<23) 

-  iVrtuJATE  AND  PRINT  THE  CONCENTRATIONS  Or 

NITRATE  AMI  HYD30X7  IONS  A  2  EACH  POINT 
IN  X  DIRECTION  AT  THIS  TINE - 


v 

C 

c 

c 


CALI  CONCEN  (C 1,C2, NX  ,N  D  ,D  ET  A,  HOOT,  DIF  1  , Y) 

■iPITE  (6, 230j  (Cl  (  J)  ,0  =  1,  NX) 

WRITE  (6, 2-j  J)  (C2  (J)  ,J=1 ,  NX) 

- CEECX  IF  THE  I'XTZSAAilON  TI.1S  AT  THIS  N  0  A  ENT  SEACnZS 

I  ti  L  NAXNIUN  TINE  TO  32  IN  1 2GB  AT  ED - 

IF  (TAU.GE.  TAUNAX)  GO  TO  1  >0 
85  !=!♦  1 

IF  (I.  E  J.  1)  Gu  TO  90 


X 

C 


-  iest  if  t;ie  co  nce  nt  ha  ti  on  g?  nitrate  ion  at  the  nth 

COLLOCATION  POINT, 12,  THE  LAST  C N 2  BEFORE  THE  SPLINi: 
POINT  IS  I  IT HI N  O-OGOI  DxNENS ION  LESS  CONCENTRATION 
UNIT  OF  THE  BOUNDARY  CONDITION.  - 


C 

c 

C 


IF  ( (  1. -Y  i  1,  )  .  IT.  0.  0001)  NO  TO  90 

-  IF  TaS  TEST  FAILS,  THE  SPLINE  POINT  *3  I  NCI.iAJ  2D,  1 2, 

NOVEL  FUPTHEP  INTO  THE  SOLJTION. - 

C  A  I*  L  «•  H  <i  N  0  Xf  (  wS ,  Di.  a ,  CO  x  r  3  ,  w  i  4  ,  CO  F5  ,  rAxiOR,DTfiU,DxAJ*"H,DIAJf'1A, 
1  J  START) 


-  EVALUATE  THE  CONCENTRATIONS  OF  NITRATs,  AND  HYDGOxY 

x  0  N  S  AT  T  H  c  CO  L  LOC  A  TI  ■-/  N  POINTS  A I  T  H  x  N  x  W  COG  A  D I N  a  2  S 
DJx  TO  THE  INCREASE  0?  THE  SPLINE  POINT.  - 


CALL  EXPAND  (7,500T, DIF  1, FACT JR, NO) 


x 

C 


INTEGRATE  EQUATIONS  (50)  AND  ;31)  USING  THE  CO  NC2N  13  ATIu 
it  .  In  x  wOL  Mj  CAT  iO  N  POxM.S  AT  THIS  TI.lx  TO  0  x  T  a  x  »*  THE 
COnCENTPATIONS  AT  THE  SUCCESSIVE  II 'E  - 


9  0  CfiLtf  D  Ir  S  J  L.  ;  N  x  H  ,  7  A  J ,  7 ,  S  A  V  x  ,  DT.i  J  ,  DT  AS  N  I  ,  0  T  A  J  N  A  ,  E  PS  ,  H  r  ,  N  A  a  , 

1  x  3  3  0  3  ,  X  F  L  A  G  ,  J  S  T  A  3  2  ,  N  A  X  D  E  3  ,  ?  -i ) 

IF  (X7LAG.  LT.  1)  GO  TO  190 
x  .  rtii.lx  =  .  JO  1  *  0  T  A  J 
JT  A  U  ,1 A  =  2  .  J  *  0 1 A  U 
.f  (if  Out  i)  jO  TO  9  3 

A>r(Y(1,1)*xx.1.«OP.Y(1,x)«sjL.  I..w3.Y(1,3).ijt.1..0.'.iv’,  ,xj  ,  j  T  «  1  . 


n  n  o  o  r  i 


lij-2 


1  .3&.3Au5('i  (1,5)-1-).uP.0.J304)  Go  TO  90 

GO  TO  93 
9  3  M1=I/IC 
M2=M  1*10 

IF  (I.  E  .12.  Oi.  TAU.  GE.TA'JMAX)  GO  TO  95 
GO  TO  d5 

- CALCULATE  THE  SURFACE  CONCENTRATIONS - 

95  XX (1) =Y  ( 1,  1) 

ITSa  N=  50  0 

CALL  1ST  STM  ( A  U  X2  ,  E  2 ,  NS  I G ,  NEE,  XX  ,1  T  ESN  ,  X  A  ,  B  A  I  OS) 

XS  (1)  =  XX  {  1  ) 

XS  (2)  =  -3.  *CUEF1*  (A  (  1  ,  1)  *Xa  (  1)  *  A  (  1,  NT)  *1. 0)  /A  (  1  ,  1)  -  A  (  1  ,  NT)  /A  (1  ,1 ) 
1  *CJcP2 

DO  150  J=1,N 

XS  (2)  =  XS  (2)  -  3.  *CCE?1  *A  (  1 ,  J* 1 )  *  {  (  1 ,  J)  /  A  (1  ,  1)  -  A  (  1 ,  J+1) 

1  *X(  1,  J*N)/A{1,  1) 

150  CONTINUE 
GO  TO  60 

190  XS  II 2  (6, 240)  DX ,  SDTA  U,  N  X,  DT  AU  ,  DTJ.ME,  DMAX,  AL,  o  5,  ZS , EPS, X FLAG, 

1  E? ,  NS IG 

NciIT E  (o,  2 **5)  BETA,  GAMA  ,00.' J1  ,CJN2 

-  FDE.1AT5  FOR  INPUT  AND  OUTPUT  STATEMENTS - 

2  30  FORSAI  l/,2SX, *  N03- •, IX, 9D 1 1.4/3 OX, 90 1  1.4/30X,  90  1  1.4/30 X, 90  11. 4 

1  /3  0X  ,9  01  1. 4/30X ,  9D  1  1  -  4/3 OX,  901  1. 4/ 30 X, ->01  1.4/30X,  *0  11.4 

2  /JOX, 9D1 1. 4/30X, 9D 1 1-4/30X, 9D1 1.4/30X,  *0  1  1.4/30 a.,  90  11.4) 

240  FORMAT  (//.  5X,  •  X  INTERVAL3  1  ,D12.4,32X,  ’  INITIAL  TAU  STEP  S-i-ZE3' 

1  ,012.4/52, ' NO  OF  POINTS  IN  X  DIP®  ',I3f2'JXr 

2  'LaST  TAU  STEP  SIZE3  '  ,  D 1 2 . 4/5  X ,  '  RE  AL  TIME  INTERVAL3  ' 

3  , 010. 4, 24 X,  ' DIF U SIGN  LENGTH3  '  , D 1 0. 4//5 X , ' AL3  ' , F 7 . 2 , JOX , 

4  ' 0  E=  • ,F7. 2,25X,' AS3  ' , F7 . 4//5X, ' EPS3  ',D10.2,27X, 

5  '  K  F  L.kG3  '  ,  I3//5X,  '  E?  3  •  ,  01  0 . 2 , 27  X,  •  NS  I G®  •  ,  I  j) 

245  FOSSA!  (//,  5X,  'REACTION  ORDER  G?  NG3-  ION  =',F7.2,25X, 

1  ’REACTION  ORDER  OF  OH-  ION  3 ' , F7.2//5X, • 1ST  CONSTANT  3 ' 

2  , j12.4,25X,  '2ND  CONSTANT  =',D12.4) 

2  50  FORMAT  (/,  25  a,  'OH-  ',2  2,901  1. 4/3  OX,  9011.  4/3CX,  9  0  1  1 .4/30  A,  -oD  11.4 

1  /3  0a,  i Oil. 4/3  OX, 9D 1 1 .4/3  OX, 901 1 . 4/30 X, 9 D 1 1. 4/ JO  X , 9 0 11 . * 

2  /30a,901 1. 4/3  OX, 901 1.4/3  OX, 9D  1  1 . 4/3 0 X, *0  1 1.4/30X, *011.4) 
310  PO  5.1  AT  ( 4 1 3 , 2F  5-  1 , 2  D 1  0.  2  ,1 5  ,  F  d.  4 ) 

3  20  FORM  AT  (4012. 5, 17,  F7.  4) 

3  30  FOSS  aT  (2F  7.  ±,  20  1  2.  4) 

34)  FORM  AT (215, 01 2. 4) 

3  5 0  FORMAT  i 1 5 X ,  1  0 F 1 0 . 6 ) 

3  60  FORMAT  (//,  D  1  .  4,  1  X  ,0  1 2 . 4 , 2  a.  ,2  0  j  0.  1  J  1 5  ,  F2  0 . 7) 

37  0  FuR  MAI  (/,2a,oF2J. 6/2  X,  or  2  0.  oj 
3  30  FORMAT  (2  X,  'CJLLCCATICN  POT.,  IS:  •) 

3  90  FORM  AT  (//,  3  a,  'TIME',  1  2  X  ,  '  C  J  a  R  £  N  T '  ,  5  X  ,  '  1 0  N  CO  CE  NT  I.  AT  lu.i  '  ) 

7  00  FORMAT  (3  OX,  10F1  0.  5/3  Ja,  1 JF  1  0.  5) 

STOP 

END 
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SUSS  our  IN  4  JC03I  (N  D,  N,  NO,  N  i  ,  AL, 3E,  DIF1  ,  3 II-  2,  JIF3  ,  SOOT) 
IMPLICIT  R  E  AL*  6  (  A  -  H,  0- 2  ) 

DIMENSION  DIF  1  (ND)  ,DIF2  (N  D)  ,  D  IF  j  i  ND)  ,  ROOT  { N  D) 


S03RGUTINE  EVALUATES  THE  COLLOCATION  POINTS  (ROOTS)  OF  TtiS 
J  A  CO  31  POLY  NO  ill  AL  OF  ORDER  N  AS  WELL  AS  THE  FISSI,  SECOND 
AND  THIRD  DERIVATIVES  OF  THE  POLYNOMIAL  AT  THE  ROOTS. 

- IN  PUT  PARAMETERS - 

IN  T  EOS  A  N  D=TtiE  DIMENSION  OF  VECTORS  Di  F 1  ,  3IF2  ,  J I F  3  ,  ROOT. 
INTEGER  N  =  iHZ  DECREE  0?  THE  JACOBI  POLYNOMIAL,  I.S.,  THE 
SUNDER  OF  INTERIOR  COLLOCATION  POINTS. 

INTEGER  NO  =  D  ECI D  £3  WHETHER  X  =  0  IS  INCLUDED  AS  A  COLLOCATION 
POIimT.  SO  MUST  3 E  SET  E^JAL  TO  1  (INCLUDING  .(  =  0) 
OR  0  (EXCLUDING  THIS  POINT). 

INTEGER  N1=AS  FOR  NO,  BUT  F02  THE  POINT  2=1. 

REAL  A  ,  3  Z=  III  E  PARAMETERS  OF  THE  J  AC  031  POLYNOMIAL. 

-  OUTPUT  PARAMETERS  - 


REAL  ARRAY  EDO T=OSE- DI ME S SI  ON AL  VECTOR  CONTAINING  ON  Ex  II 
THE  N*NQ  *  N1  I  EROS  OF  THE  POLYNOMIAL  USED 
IN  THE  COLLOCATION  ROUTINE. 

REAL  ARRAY  =  T  HR  E  E  ONE-DIMENSIONAL  VECTORS  CONTAINING 

DIF1.DIF2, DIFs  ON  EXIT  THE  FIRST,  SECOND,  AND  THIRD 

DERIVATIVES  OF  THE  POLYNOMIAL  AT  THE  LEGOS. 


;  ********  *  *  *  *  **£ 


C  -  FIRST  x VALUATION  0?  COEFFICIENTS  IN  RECURSION  FORMULAS  - 

C  - RECUAjIGN  COEFFICIENTS  ARE  STORED  IN  DIF1  AND  DIF2  - 

C 

A  3= A L* BE 
AD=UE-AL 
A?=3  E*  AL 

Di.F1  (1)  =  (AD/ (AB*2)  *1 )  /2 
DIF2  (T) =0 

IF  (N  .LT.  2)  GO  TO  15 
DO  12  1=2, N 
*  1=1-1 
~=A3  *2-*-1 

DIF  1  (I)  =  (AD  *A  D/2/  (2*2)  ♦  1)  /2 
IF  (I.  NE.  2)  GO  TO  1  0 
DIF2 (I)  =  iAj*rt?*L1) /2/2/ (2*1 ) 
ju  TO  12 
2  =  L*  L 

Y=  21  * ( A3*  *>  1  ) 


10 


c 

c 


10-4. 


t=i*  (A  P  ♦  I ) 

0x72  (I)  =  i'/~/ ii-1) 

12  CO 21  IN J  E 

- - -  DOT  27 I  2A  TI  0  ‘1  3?  XZitTON  .'IZItOu  WITH  S'JSPSFS SI  02 

OF  PLaYIOUSLY  OFT  23  .11  2  £  J  auOTS  - 

15  x=o 

Du  35  1=1,2 
20  .<D=0 
XN=1 
:<D1=0 
X  S 1  =  0 

t/O  2  x  J  =  1  ,  *■» 

X?=  (0171  (J;  -a)  *XN-CIF2  (J)  *X„ 

X?1  =  (DIF1  (J)  -X)  *XN  1-DIF2(J)  *XD1-XN 
XD=X2 

ad  i=x:m 

X  S  =  X  P 

22  x:n=x?J 

2C=  1 

2  =  x2/X2  1 

IF  (I .  £*.  1)00  ro  3  0 
DO  2  5  J  =  2 , 1 

25  2C=LC-Z/ (X- SOOT  (J-  1)  ) 

30  2=2/20 
X=  X-  - 

I F (DAB S  („)  .  01.  1.D-C9J0C  TO  20 
FOOT  (I) = X 
X  =  X».000  1 
35  continuz 

- ADj  2  V  £  2 1  u  A  L  COLLOCATION  POINTS  AT  X=J  OS  a=1 - 

2T  =  N  +N0  42  1 

a*  ( 2  0 •  c* o  •  '  < )  ^  0  TO  4  5 
DO  40  1=1,2 
J  =  .i  *  1“  I 

UO  SOOT  (J  4  1  )  =  2  OuT  (J) 
aOOI  (  1)  =0 

45  IF  (:i  1.  £*.  1  i  20 or  (NT  )  =  1 

- 20.»  £7  A  L  J  AT  E  DEtiIVATI7ZJ  OF  POL  t  20.11  AL - 

Do  5  0  I = 1 , 2  T 
I  =  r.JoI  (I; 

DIF  1  (1/  =  1 
0 IF2  (i>  =j 
JIF3  (I) =o 
DO  o  o  J  =  1  ,  .i  1 

a-  (  j  « i  ^  }  j  w  i  0  3  0 

£-a~ aos  r  (j) 

DIF3  (I)  =1*01 1-  3  (I)  4  3=0172(1, 

0112  (I,  =  L  *  Ji:  2  (I)  «■  2*  DIF  1  (I) 


t 


( i  n  <  i 


14 


3  IF  1  (I)  =x  *D  lx’  1(1) 
50  ccnti.nje 
fztu  a:i 
end 


suuaouriwz  jfopr  (n d,  n,  no,  n  i  ,i,  id,  dif i ,  oifi,  diz 3 , rogi,  v z. c r ) 

1.1  PL  10  IT  RZAL*d  (A-H,0-  Z) 

DIMENSION  DIF1  (HD)  ,  DIF2  (HD)  ,DIr  J  (HD)  ,2(001  (HD)  ,  '/2CT  (HD) 


:***«#*****#: 


c  sajkouriHE  evaluates  discretization  coefficient  maijices  c 

<:  AND  GA  USE  J.  AN  wUADRATURS  WEI  Gil  ID,  N0.1ALI  LED  TO  S  J  .1  1-  C 

c  c 

C  - INPUT  PARAMETERS  —  c 

c  c 

C  INTEGER  J,  JO,  N1=AS  IN  SUBROUTINE  JCCDI.  C 

C  INTEGER  I=THz  INDEX  CF  THE  NODS  FOj  WHICH  THE  WEIGHTS  ASS  z 

TO  ot  CALCULATED.  C 

INTEGER  ID=1N DICATOR.  ID=1  GIVES  THE  WEIGHTS  FO 3  Di/DX,  C 

ID  =  2  GIVES  THE  WEIGHTS  FOE  D2T/JX2,  AND  ID=J  <- 

GIVES  THE  GAUSSIAN  WEIGHTS.  THE  VALUE  OF  I  IS  C 

C  la  RELEVANT  IN  THE  LAST  CASE.  C 

c  seal  aprat  root,  c 

C  Jl F 1 , D IF2 , 0  IF  3  =TH2  ONE- DI ME NSIONAL  VECTORS  COG PUT  ED  IN  C 

C  SUBROUTINE  JCOSI.  C 

C  c 

C  -  OUTPUT  PARAMETERS  -  C 


C  REAL  A  a  3  AT  VECT=THS  COMPUTED  VECTOR  OF  « EIGHTS.  C 

c  c 

NT=N +N0+N 1 

IF  (I  D.  e;.  3)  Go  TO  1  0 

DO  20  J  =  1  , J  i 

IF(J.NE.  I)  GO  TO  21 

IF  (I  0.  nZ.  1 )  Gu  TO  5 

VECT  (I)  *  Di?  2  ( I)  /D I  FI  (I)  /2 

Gu  TO  20 

5  VECT  (I)  3  Jxf  J  ( I)  /D I  r  1  (I )  / 3 
GO  4  o  20 


21 

T- AO OT  ( I ) 

-LOOT  (J) 

J  ECT  (J)  =D 

IF  1  (I)  /DIM  (U ) /£ 

i;  1 1  D.  z  *  • 

2)  VZCT(J)  =VZC  T  (J)  *  (DIF2  (I)  /JIF1  kI 

r  > 

CwNliNUZ 

J  J  4  0  3  0 

10 

D'.y  2  5  J  =  1 

,<  r 

::  =  Fo^T  ijj 
A  X  =  X  *  l  1  -  a  ) 

i  F  (  N  0  .  Z  .  .  J)  A  a  —  \  X  /  X  /  X 
*  x  (.1  1 .  x*  a  *  0 J  *i X  —  A  a /  (1—  X)/(1—  a) 


Uo 


VZCI  (J)  =  A  X  /  J I  F  1  (J)  +*  2 
25  l  =  Y  ♦  /  E  C I  [  J  ) 

CO  o  0  J  =  1 ,  N  T 
60  VECI  (J)  =  Vr.CT  (J) /Y 

50  return 

END 


s  ubrootine  i.its?  (n  c,  nt ,  x,  soot,  di?i , xintp) 

IMPLICIT  R2aL*8  (A-H,0-Z) 

DIMENSION  ROoT  (ND)  ,DIF1  (ND)  ,XINT?(ND) 


C  *****  1*  I 


(«*****«* 


c 

C  SUBROUTINE  EVALUATES  THE  LAGS  ANSI  AN  I  NT 


******#************»****(_■ 

c 

SR  PoLA  710  N  OOLFF IOIZNTS 


C  - INPUT  PAR  AHSTE3S -  C 

C  L 

C  INTEGER  N  I  =  I  tt  Z  TOTAL  N'JHjER  OF  COLLOCATION  POINTS  (  =  N +:.0  ♦  1 )  C 

C  FO  a  WHICH  THE  VALUE  0?  THE  DEPENDENT  VARIABLE  Y  C 

C  IS  KNOWN.  C 

C  REAL  X=TtiZ  A3SCISSA  X  WHERE  Y  (X)  IS  DESIRED.  C 

o  PEAL  ARC  AY  C 

C  SOOT,DIF1  =COLLOCATICN  POINTS  AND  DESIVhTIVES  CF  NODE  C 

C  P  J  u'i  NO;".  I  A  L,  DERIVED  IN  GU3RO  UTINZ  JC'j3I.  C 

0  C 

C  -  OUTPUT  PARAMETERS  -  C 

0  C 

C  SEAL  ARRAY  XI  NT  ?=T  H  E  VECTOR  OF  INTERPOLATION  WEIGHTS.  C 

r-  C 


*******  i 


:**♦¥** 


PCL=  1 

DO  5  I  =  1  ,  N  T 
Y= X- ROOT  ll) 

XINTP  (I)  =0 

IF  (i  .  2*.  0.  00)  XINTP  (I)  =  1 
5  ?GL=  POL  *  Y 

IF  (POL.  S*.  0. 00)  GO  TO  10 
DO  o  1=1, NT 

o  XINTP  (I)  =POL/DIF  1  (I)/  (X-30OI  1 1)  ) 
10  RETURN 
END 


S  J  J  ti  u  U  T  i  .i 

1 

IMPLICIT 

Jl.-iENSIDN 

3 

0  x  X  o  ft  S  I J  N 


E  Oxfcj'lo  (*^,  A  /  . ,  d  ,  >i>.  x  S  ,  «1  ■!  .C  ,  wT  j  f  .i  r  ,  i 

JS TA  R  T ,  U  A  X D E  R  ,  P  .1 ) 

?  u  ilL  *  J  (A™3,  ^  ~  w  ) 

Y  (0,20)  ,  Y.l  AX  (  20  )  ,  SAVE(2A,xO)  ,  ERR  OS  i  2  0  ) 
A  to  i  ,P'£STST(7,2,j) 

«.  UJ)  ,?  JL(20,20)  ,aUJ) 


-1ft 


ft , 


kl  J  <J  D 


,?»  k  -•  0  J  j  , 


«T  L  U  ; 


«*«*  ***»*£i***4*>a<i*w#**i»*^*^ 


o  n 


147 


\ 


\ 

i 


i 


c 


jilJS'JJ  Cl HE  xNlEGRAT 
ORDER  EQUATIONS  AS 
ONE  STEP  0?  LENGTH 
DECREASED  .tliulN  TH 
AS  LaaGE  A  STEP  AS 
STEP  SHiiDR  #  a  IC  !1  IS 
EACH  C OH PONS NT  G?  7 


ES  A  SET  0?  .<  ORDINARY  DIFFERENTIAL  FIRST 
DESCSI  JED  ID  EQUATIONS  (50)  AND  (51)  OVEc 
H  AT  EACH  CA^L.  a  HAY  HE  INCREASED  OR 
E  a  A  USE  1U.1  TO  HMAX  IN  ORDER  TO  ACHIEVE 
POSSIBLE  *  SHE  NUT  COMMITTING  A  SINGES 
LARGER  THAN  EPS  IN  THE  L~ 2  NOSH,  MHxLE 
HE  ERROR  IS  DIVIDED  BY  THE  C0A?J:,Z.,Z3  of 


YM  AX. 


TH  S  PROG  RAH  CALLS  FOSS  SUBROUTINES  NAMED 
DIFFUN ( T  ,  Y  ,  D  T ) 

:i  A  T  i  N  v  (Pm,n,m,  j) 

SOLVE  ( it  N  ,  c'U  L ,  3 ,  X) 

P2DERV  ( T  ,  Y  ,  ?  .i ,  .1 ) 

THE  FIRST, Dir  FUN, EVALUATES  THE  DERIVATIVES  OF  CHS  DEPENDENT 
V  ARIAS  EES  STORED  IN  Y(1,I)  FOR  1=1  TO  N ,  AND  STORES  T  ;i£ 
DERIVATIVES  IN  THE  ARRAY  CY.  THE  SECOND  IS  CALLED  ONLY  IF 
Hr  IS  SET  ID  1  OS  2  FOR  STIFF  HSI.ioDS.  IT  PEHFOR  IS  A  FlaST 
STAGE  OF  GAUSSIAN  EL  IM I M  ATI  ON  ROUTINE  OF  TiiE  N  3Y  N  MATRIX 
STORED  IN  THE  ARRAY  PM  (M , M)  .  SOLVE  PERFORMS  THE  jACK 
SIJSSTITUilON  PROCESS  0?  THE  GAUSSIAN  ELIMINATION  RoJTlN-. 

P  x.  ij  Cr  .(  V  x  j  .-j  S  o  0  0  N  u  Y  xr  Hr  I  j  1  ,  A. ID  a  .  O  a  *—  S  x  a  £>  C  /i  r  ^  in. 
DERIVATIVES  Dr  THE  DIFFERENTIAL  EQUATIONS  PROVIDED  DY  DIFFUN. 

-  PARAMETERS  — - 


N=  THE  NUMBER  OF  FIRST  OR D ES  OIFFERr-N  TIAL  EQUATIONS. 
T=riiE  INDEPENDENT  VARIA3LS. 


c 

Y= 

A  d 

BY 

N  ARRAY  CONTAIN 

ING  THE  DEPENDENT  V 

ARIA3LES  And 

c 

W 

the: 

h  1>Z 

ALED  DERIVATIVE 

S.  THE  DEPENDENT  V  A  SI  A3  LEE  Ahc- 

c 

z 

3  a  o  £1 

ED  I 

I  (1,1)  ,  l=i  T 

C  N. 

V- 

* 

wi  A 

V  E=J 

BLO 

CX  0?  AT  L  £A  ST 

12  *N  FLOATING  POINT 

Z  w  A I  i.  C  Z  £  J  j  u  j 

c 

c 

b 

Y  TH 

a  SUBROUTINE. 

c 

H= 

THE 

3TZ  ? 

Si  Gw  TO  3  E  ATT 

SMPTED  ON  THE  NEXT 

STEP. 

w 

C- 

HM 

IN=T 

HE  MINIMUM  STS?  SIE 

E  THAT  MILL  SE  USED 

F 03  THE 

c. 

1 

»i  r  z  j 

RATION. 

C 

:ih 

a  a—  r 

HE  M 

A  X i M  U  M  SIZE  TO 

MHICH  THE  STEP  MILL 

3  w 

c 

S?S  =  Tii 

Z  £  it  n  Z  a  *  Z  Z  T  CCNSTA 

NT- 

C 

J 

=  THE 

;isr 

HOD  INDICATOR. 

Vw 

0 

AN  ADAMS  PREDICTOR  CORRECTOR  IS 

USED. 

c 

Z 

1 

A  MULTI-STEP 

METHOD  SUITABLE  FOa 

STIFF  S  Y  S  T  E  G  S 

Z 

IS  USED. 

c 

J 

Z 

THE  SAME  AS  C 

ASE  1,  EXCEPT  THAT 

T  H I  S  S  U  3  K  g  U  1 1  N  - 

z 

■* 

COMPUTES  THE 

PARTIAL  DERIVATIVES 

J  :  xi  KUaiw.u 

z 

£ 

DIFFERENCING 

OF  THE  DERIVATIVES. 

c 

z 

7  •/ 

A  X  =  A 

.1  A  a  a  <W  u  r  ?«  I*  CC  AT  I 

ON S  Mu I CM  CONTAINS 

4.  .lw  ■!  .i  ■!  l  kl  J  |V|  V>  fc 

z 

j 

5 

ACH 

Y  SEEN  SO  FA?. 

z 

EZ 

ROa= 

AN  A 

PRAY  0?  M  ELSMZ 

.<?  3  *»G  IJT  j.  .1 Z 

X  '.W  X  4  .4x14 

c 

u  A  Z 

.IF  St-  ERROR  IN  E 

ACM  CQ.lPONENT. 

c 

V  ^ 

.1  xl“ 

A  CD 

M  P x. E  T ION  CODE. 

^  xi  £,  A  Z  £  l  ii  w  *  -G  Z  ^  o  juoLwJ^r  J  Lj  + 

z 

J  5 

r  \  a  r 

=  AN 

laps:  r;  o  i  ca  to  3 

« 

z 

z 

. .  it 

X  D  ER 

=  T  d  E 

ijAaImj-i  d e?  : v  \ 

::v z  :.ur  z .  i  z  j  z  z  3 z 

u  3  -  J  • 

V. 

c 

PM 

*  A  b 

uO  wA 

OF  a  T  LEAST  N  *  *2  rLO.iTINO  pu  I  N  T  LOCA  T  .0  NS. 

z 

1 


*«***«**4.**««****«*.»*#**(2 
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********«*****?>  V«**4*9  *  *  **  *  *  *  *  * 


DATA  P  E  H  x  J  T  /2.0,4.5,7.33j,10.42,13.7,17.15,1.J, 

1  2. 0,1 2.  3, 24. 0,37. 69,  53. 33, 70. 0  6, a  7. 9  7, 

2  3.  0, o. 0, 9.  1  o7,  12.5 , 1 5. 9d  ,  1. 0, 1 . 0, 

3  12.  0,^4. 0, 3  7. d9,  53.33,  70.06,3  7.97,  1 .0, 

4  1 . ,  1.  , 0.5, 0.1667, 0.04133, 0.0 0d2o7, 1.0, 

5  1.0, 1.0, 2. 0,1.0, .3157, .07407, .0130/ 

DATA  A  (2)/-  1.0/ 

I H  ET  =  1 
KFLAG  =  1 

I?  I JSTASr. LS. 0)  DC  TO  140 


1  oo 

DO 

1  10 

1=1,3 

DO  1 

10  J=  1,  K 

1  10 

3A  V  E  (  J  ,  I )  - 

i 

(0,1) 

:iULD  - 

a  0  2  « 

IF 

(a. 

EQ. HOLD) 

GO 

TO  1  30 

1  20 

£<  acj  ft 

=  11/  3  CL  D 

13 

ET  1 

=  1 

GO 

TO 

750 

1  30 

ilQOL  D 

=  N  w 

T  Oi.D  =  I 

k\  A  C  U  *1  —  i  •  0 

Ir  (J3TAbT.Or.Oj  GC  TO  250 
GO  TO  170 

140  IF  (JSTA.iT.  E  *  .  -  1 )  GC  TO  1o0 
-«0  =1 
:i3  =  n 

;n  =  n  •*  i  o 

N2  =  N 1  ♦  1 
;;  4  =  N  **  2 

:«5  =  ni  ♦ 

:io  =  ii  5  *  i 

CALL  OIF?  J  (T,Y,3A7E(N2,1)) 
DO  150  I  =  1,3 

1  50  i(2,I)  =  JA7Z(N1+I,1)  *.i 

’.IN  2*  =  J 
K  =  2 
30  TO  100 

160  IF  (  N£.  2'j.  il  vG  LD)  JSTABT  =  1 
T  =  TOLD 
NC  =  N^Oi.0 


K  - 

♦  i 

GO 

TO 

120 

1  70 

IF 

(IF 

•  £»  ^  • 

■j  0 

TO 

1  30 

1  F 

(Jv 

-  g  r  • 

=>  j 

.j  G 

TO 

1  GO 

GO 

TO 

(22  1 

■) 

/  «- 

22, 

223 

,224  , 

22  3,  22  o )  , * 

1  30 

L? 

(  by 

•  J  A  • 

7) 

GC 

TO 

1  GO 

TO 

(2  11 

»  «• 

12, 

213 

,214, 

215,216,217) 4 

1  30 

Kr  L 

A  a 

=  “ 

f.  57 

•ja:i 

2  1  1 

All)  = 

-  1. 

G 

GO 

:u 

2  30 

2  12 

a  r- 

)  = 

5  j 

J00 

00000 
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A  iJ|_ 

=  -  •} .  5  o  j  0  00  0  0  0 

JO  k 

u  1 3  J 

u 

13 

A  (1) 

=  -J.  -t  loo  6 06060006 6o 7 

A  l  3) 

=  -0. 750000000 

A  ('*) 

-  -j.  1  o  06  6  o o  6  6b  66  fa  o7 

jO  T 

0  Z  30 

2  14 

A  (I) 

=  -0-375000000 

A  ( 3) 

=  -  0- 9  1  06 b  co  6  6 666  6  66  7 

A  (4) 

=  -0.  3  33332  33  33  33333  3 

A  (5) 

=  -0.04l6666666  66  6  66  6o-0  7 

30  10  230 

2  15 

A  ( 1. 

=  -0. j4oo  1  1  1  1  1 1  1 1  1  1 1  1 

A  iJ) 

=  -  1 .  0  4  1  66  6  o  6  66  66  7 

A  (-) 

=  -0. 4  Jo  11  1  1  1  1 1  11  1  1 11 

A  p) 

=  -u.  10 4166 606666660 7 

A  (o) 

=  -O-OOo  3333  3  3333  333  3  3  J 

jO  r 

0  2  3./ 

1 

1  6 

A  ( 

=  -  0. 32* 36  1  1  1  1 1  1 1  1  1 1  1 

A  13) 

=  -1.1  -416066660060  66  7 

A  ( 4) 

=  -J.o2o00000C 

A  (5) 

=  -0. 177 )3 33 333333 33 33 

A  (o) 

=  - J. Jz 50000000 

\0\ 

=  0.  00  1  3  3  33  6  3  S  d  d  3  3  36  60 

Go  r 

0  230 

“> 

<- 

17 

All) 

=  -u. 3  1o59 16  3  12  16931  2 

A  l  3) 

=  -  1.  ZZ50J000C0 

A  (4) 

=  -0.  7  5135135  1351  351 J 

A  (5) 

=  -O.265206333333  33333 

A  l») 

=  -0.0  436  1  1  1  1  1  1  1  1  1  1111 

A  (7) 

=  -  0.  0  04  06  1  1  1  11  11  11  1  1  1  1 1  1 

A  (o) 

=  -0. 000  19941  26934  12o9d4 

o  0  r  J  230 

2 

2  1 

A  01 

=  -I.OOdOOOOOuO 

4. 

0  230 

222 

A  (1) 

=  -O.00006 6 666600606 7 

A  i  3) 

=  -0.  3j  03  333  3  33  33  333  3 

r*  •** 

J  * 

U  23  0 

■) 

A. 

23 

A  (  1) 

=  -j.  5454545454 5**5-04 5 455 

A  (3) 

=  Ail) 

A  (4) 

=  -O.u0090909G9090  90 

3o  r 

0  230 

T 

24 

A  (  1) 

=  - J.40OOOOOOG 

A  ( J ) 

=  —  0  .  /  00  0  00  0  0  0 

A  (  4 } 

=  -0. 2o0000000 

A  P) 

=  -u. J20J009JC 

JO 

0  z3  j 

■) 

z5 

A  k 1 ) 

=  -0.4 37 95020 43 79562 

A  (3) 

=  -  ).  jZl1b7jd32  1lo73J 

A  i  4  J 

=  -j.  3 luZ 1 5Q7d102 1 34j 

A  (  5) 

=  -  0 .  J  5-»  7  44  5 2  55  < 7  4 45  2  c 

— 

-3  •) 

4  .  / Jj_n^  i  6  3  5J  36  4  06  j  3  J  4 

JO 

10  Z  JU 

1 

2  c 

AO) 

-  -j.. do  1 6  32  6  53  Jc  1225 

A  Pi 

-  -u.  .#  zO  0  3  4  5  2  Co  34  9  Z  0  0 
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A  (4)  =  -0 .  4  1  o  6 6 06  5  66 oo  o uti  o  7 

A  (5)  =  -0.  03*2047619047619 

A  (t)  =  -0.  0  1  193476  1904  761  j 

A  1 7 )  =  - 0 .  0  00  566  3 9  34  240  Jo 2  3  2 
2  30  X  =  NQ *  1 
I DuU B  =  K 
Ml'Y?  =  (4  -  ;i7)/2 
E3Q2  =  .  5/f,‘  LO  AT  (N  Q  ♦  1) 

BNOi  =  .  5/F  L J  A  T { N  0  ♦  2) 

ZNW1  =  .  5/T  LOAT  (:JQ) 

PEPoh  -  ZP3 

EU2  =  ( ?  Z  3  2  j  i  ( 3Q , B  TY  ? ,  2  )  *  ?  E  2S  d )  *  *  2 
3  =  (PSaiST  (0C  ,KTYP,  1)  *?z?Sd)  **2 
2DAN  =  ( 2  Ea 18  T  (0  0#  .IT  Y? ,  3)  *2  EPS  H )  *  *2 
IF  t2D«N.2Q.O)  GO  TC  780 
3N  0  =  EP  2*i 30  3/DFLCA  I  (M) 

2  40  I  ii  Z'J  AL  =  111 

go  ro  (250  ,odOj  ,  ihet 
2  50  r  =  i  tl 

DO  2  bO  J  =  2 ,  X 
DO  260  J1  =  J,K 

J  2  =  X  —  J  1  *  J  -  1 
DO  2c 0  I  =  1  ,  '3 

2  60  i  ( J  2  #  I )  =  Y  ( J  2  ,1 )  *  Y  ( J  2  *  1  , 1 ) 

Do  2  70  x  =  1,3 
270  EaaOa(L)  =  0.0 
DO  430  L  =  1,3 

CALL  DIFFU3  ( T  ,  Y  ,S  A  VS  {  32 , 1 )  ) 

IF  ( I'«  2VAL.  LT.  1)  GO  TO  350 
IF  {»?.£*'.  2)  00  TO  310 
CALL  2  E  D  2  3  V  (T ,  Y  ,  Fi»  ,  3  3  ) 

A  =  A  (  1 )  *  A 
DO  280  I  =  1 ,  M 4 
280  ?3  (I)  =22  (I)  *a 

290  00  300  I  =  1,  N 

3  00  p  w  ( i » ( a  3  ♦  i )  -  a  3 )  =  i.o  +  ?*  (i  *  (ai*  i ) -;ii) 

I  4  EV  AL  =  -  1 

CA  ..  (1ATI3  V  3  ,  3  3,  J1  ,  PUL) 
i.F  (J1.GI.0)  GO  TO  3  50 
oO  TO  440 

310  DO  320  1  =  1,  N 

320  S  A  V  -  ( 9 , 1 )  =  Y(1,I) 

do  340  j  =  i ,  :j 

a  =  22  0*  D.IAa  1  (E  ?S  ,  DA  35  (SA  V£  (9  ,  J)  )  ) 


Y ( 1 , J)  =  Y  (1,  J)  ♦  ? 
J  =  A  (  1 )  *;i/n 
v.  A  L ..  3  .  .  r  J  .1  (  4  ,  {  ,  j  A  / 
DO  3  j  0  I  =  1 ,  a 

z  ( a  6 ,  i )  ) 

3  30 

2.  (I*  (J-  1)  *>13)  = 

(5  A 7  £  ( :i  5  ♦  1 ,  1 )  -  0  A  V  2  (a  1  *  1 ,  1 ) 

/  *■> 

3  40 

Y  ( I,  0)  =  2A  7  2  (  S,  J) 

GO  To  2  )<) 

3  50 

IT  1 2 2. a  2.0)  GO  ro 
DO  3  o  0  I  =  1  ,  3 

3  70 

360 

3AV2(*,I)  =  Y  ( 2 ,  I 

)  “  Gu  /  i  (Ii  i  *1  ,  1 )  *  :i 
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CO 

X  D  4 

1  J 

3  70 

DO 

3  60 

I  = 

1  ,  N 

3  30 

CAVE 

( N  5  ♦  I 

,  1) 

=  Y  (2,  I)  -  SAVE  (N  UI,  1)  * H 

call  s 

Du  VE 

(•< 

,?  J  L  ,SAVS  (N  34  1  ,  1  )  ,2) 

DO  399 

0=1  , 

N 

3  99 

SAVE  (9 

,  J)  =.c 

(0) 

4  10 

NT 

=  N 

DD  4  20  I  =  1,'l 

UM)  =  t  (1,I)»A(1)*3A/S(9,I) 

Y(2,I)  =  1(2,1)  "  SAVZ(4,I) 

33*02(1)  =  S3  ?  CP  (I )  ♦  CAVE  (3,1) 

IF  (DAdS  (3A  73  (9,  I)  )  .  L3.  (3ND*Y.1A.C  (I)  )  )  t. T  =  N  I-  1 
4  20  CONTINUE 

IF  (DX.L2.0)  CO  TO  4*0 
4  30  CONTINUE 
440  T  =  TOLD 

IF  t  (H.  L2.  (  HdI'J«  1. 0000  1  ) )  .  A  ND.  {  (I  WEVAL-MTY  P)  .  LI. -1)  )  uo  T 
IF  (  (SF.  tw.  0)  .03.  (I«  ZVAL.  NS.  0)  )  3ACUM  =  a  ACU  1*  .  2  5D0 
IWEVikL  =  he 
I3ET1  =  2 
CO  TO  7‘j0 
^  60  KFLAC  =  -3 
4  70  DO  4  d  0  x  =  1  , 

DO  4  d  0  0  =  1  ,  K 
4  60  Y  t  J  ,  j.  }  =  CAVE  (J,I) 

F.  =  HOLD 
JQ  -  NOOuD 
JSTA3T  =  NO 
RETURN 

4  90  D  =  0-  J 

DO  5  00  I  =  1 ,  .! 

5  00  D  =  J  *■  (Z3R03  (I J/  YMAX  (I)  ) **2 

I  *  c. V  AL  =  0 

IF  (D.  CT.  E)  CO  TO  540 
IF  (rC.LT.J)  CO  TO  520 
DO  510  J  =  3, K 
DO  5  10  1  =  1 ,  N 

5  10  Y  { J  ,  I )  =  Y  ( J  ,  I )  ♦  A  (J)  *E53U3  (I) 

5  20  KFLAC  =  f  1 
HtJE'n  =  II 

IF  ( iDO U B . u E. 1 )  CO  TO  550 
IDOU3  =  ID0J3  -  1 
IF  {  IDOU6.  C  X.  1)  30  TO  700 
00  5  30  I  =  1 ,  :i 
530  3*  7E  t  10,  I)  =  333CE  (I  ) 

CO  TO  / 0 J 

540  XFLA  j  =  iC*  Ui\u  ~  2 

IF  ( H.  u2.  (ii  AIN  *  1 . 0  000  1 )  )  Co  TO  740 
T  =  TOLD 

i-*  (KF uAu*  i#i«  *5)  jC  TO  720 
5  50  ?P 2  =  I D/-;  **2-;;2*  1.  2 
Pa  3  =  1.;*J 

.  *  (  i  a  2  *  4  u  .  H  «  a  D  &  E )  .Oa.  (  Ka  ua  j*  utii  •  ~  1 )  J  TO  3  7  J 

D  =  0.  J 


0  4  o0 


1 
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DO  5  oQ  1=  1  ,  N 

560  D  =  D  *  (i2550a{I)  -  SA  72  (1  0,  I)  ) /YliAX  (I)  )  **  2 
PSJ  =  (D/ZJ?)  «*ENQ3*1.  4 
570  Pal  =  1.2*20 

IF  (NQ.  L£.  1  )  GO  TO  530 
D  =  0.0 

DO  5  80  I  =  1 ,  a 

5  80  D  =  D  *  (  7 (K,I)  /YJ1AX  (1)  )  *  *2 
?3 1  =  ( D  /  E  J  «i  >< )  **EMQ1*1.  3 

5  90  CONI  IN  JE 

IF  (PP2.LE.P33)  50  TO  6  30 
IF  (P3J.LT.  Pa  1)  GO  TO  660 

6  00  3  =  1. 0/.\:1Aa1  (?R  1  ,  1.  E-4) 

N2.»Q  =  NO"  1 
6  10  I  DOL'D  =  10 

I?  (  (KFLAa.  2*.  1)  .AND.  (R.LT.  (1.  1)  j  )  GO  TO  700 
IF  { N E n * . L E . N i )  GO  TO  630 
DO  o 20  I  =  1,N 

6  20  7  (N2*2*1  ,  I)  =  28  30  8  ( I)  *  A  (  X)  /DFLG  A  T  (X) 

o  30  K  -  N  E  ’<i  ^  *  1 

IF  (KFuAG.  L*.  1)  GO  TO  670 
aACJd  =  aA2J>l*S 
I  a  21  -  3 
GO  TO  7^0 

640  1?  {  N E «  j .  .1  v )  GO  TO  2  50 

M  Q  =  >12  W  j 
GO  TO  170 

650  IF  ( PK2.G1.  PL  1)  GO  TO  600 
NENQ  =  Nw 

P.  =  1.  J/AMAX1  (?R  2  ,  1.  E-4) 

GO  TO  olO 

6  60  R  =  1.  J/Ail  A  X  1  (Pit  J  ,  1.  E-  4  ) 

NEWQ  =  N  ^  ♦  1 

GO  TO  610 
670  13  EX  =  2 

R  =  DM  IN  1  (5  ,  i(8  AX/D  A3  2  ( H) ) 

ii  =  u«y 

LNErt  =  H 

IF  (  N<*.  2Q.  N  £  *  2)  G  O  TO  6  30 
HQ  -  NSNg 
GO  TO  1  7  J 
680  5  1  =  1.0 

DO  6  30  J  =  2,  X 
21  =  31*6 
DO  6  >0  I  =  1  , 

6  90  7(3,1)  =  Y  (0,1)  *31 

I  JOUR  =  K 

700  DO  7  10  I  =  J , N 

710  l.!AX(I)  =  J  MAX  1  (  Yfl  AX  ( I)  ,  D  Ac.  J  (  7  ( 1 ,  I)  ;  j 
J  3TA  3T  =  2  * 
i  2T’J  NN 

720  17  { N Q . E ^ .  1 ]  GO  TO  7  oO 

C  Ac.  DiFFJN  (T  ,  7, SAVE  (N2,  1 )  ) 
a  =  H/riOcJ 
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do  7  jo  :  =  i ,  :i 

7  (  1  ,  1)  =  SAVE  (1  ,  I) 

SA7E(2,I)  =  :iGLD*5AYE  (NCI,  1) 

7  30  i  [2,1)  =  J«Vi(2  ,  I)  *3 
0  =  1 
K  F  LA  G  =  1 
00  I  0  1 7  0 
7  40  KFLAO  =  -1 
KNZrf  =  fi 
J  START  =  N2 
SETH  DN 

750  :<A  Cod  =  DHAXipABS  (H  HI  U/H  OO)  ,  3AC'vM) 
EACH  0  =  OH  I  A  1  (3  AC'J  M, DA  3S (  371  AX/HoL  J)  ) 
Si*  1-0 
DO  7  60  J  =  2,  '< 
n  1  =  H  1  *a  AO U H 
DO  7  o  0  1  =  1  ,  N 

760  7(J,D  =  SAVE  (J,I)  *51 

ii  =  holo^aacuh 
60  7  70  I  =  1  ,  N’ 

7  70  1  (  1  ,  1)  =  SAVE  (1,1) 

I  jjO  0  0  =  i\ 

00  TO  {100,250,640), ID ET1 
7  30  KFLAG  -  -4 
GO  TO  47o 
J 


SJBBOUTIHE  DIFF'JS (T, YY  ,0YY) 

IMPLICIT  -iZAL*6  (A-H,  j-Z) 

DXG&SSUM  0  YY  (23)  ,  YY  (3 ,20)  ,  A  (20,20)  ,3  (20,20)  ,C  (20)  ,.(2(1)  ,W  A  (3) 
1  ^  A3  (20) 


■  •  * **  *  * i 


o  sjl3wujti.j2  called  3y  cifsob  evaluates  the  derivatives  of 

C  D  2  ?  Eli  D  23  I  VARIABLES  STORED  1.1  Y  L  (  1  , 1 )  ,  1  =  1  ,  2  *  A  ,  .» :i  D  STORED 

C  THE  DERIVATIVES  ill  THE  ARRAY  DYY.  I  IS  THE  Hi  LED  EN  DEli  I 

v.  /  A  i\  X  A  3  u  &  • 


C  **  3 


i  ********  i 


‘C 


L  a  T  3  a  N  A  I*  A  J  X  1 

^N  /  ul  /CuilFI  /  CO  ZF5 
CJVA  0:1/12/3  0  2?  2 
v  3  3  C  !<  /  l.  3  /  .•  t  k*  T 
CO  *1.1 ON  /L 4 /A  ,  o 
cj.u  j;«/l7/:.31-3,  ;i2  3,Z? 

:i=j  »:i 

4*  J  3  0  I s  i  ,  I  I 

«*a  (n  s  i  *  l  ^  i  i) 

30  v.  j  ti  i  ^  N  J  w 

— —  ■“  j  i  *  H  i  /A  w  Jwo  \j  3  x.i-i  3  J  «?  Ao  -  ^a—)) 


15^ 

ii(l)=YTll,  1) 
i  r  es  n=  5  ou 

ll  *«U  4j  *J  4  ~J  X  »  i  ^  ,1  J  .C  1  /  H  r*  ,  l  W  |  I*  A  ,  .1  u  ,  X  3  ) 

CT=-  J.  *CjEE  1*  (A  (1  ,  1)  *A.i  (1  )  +  A  1 1  (j;)  *1.  ,  /A  (1  ,  1)  -  A  (  1,  „2)  /«  1 1 ,  1;  «o  0 
DO  4  3  i  -  1  ^  d 

C  (I  j  =  3(1+1,  1 )  *CT  +  5  (1+  1  ,NI)  *C  JaF  l 
45  CONTINUE 

—  STORa  IHZ  DEFI7AT IV  ZJ  IN  ANPiii  D'l'l  — - 

30  1  20  J  =  1  , 

DTi  (J)  =  j  (0  +  1  ,  1  )  *.{ .(  ( 1)  +  3  (J  +  1  ,  NT)  *1. 

3iY  (J  +  N)  =0  iJ) 

33  110  K=  1  ,  :j 

Dli  ( J  )  =  3  1 1  ( J )  +3  ( J  +  1  ,  X  +  1 )  *'i  i  (  1  ,  K  ) 

3  Vi  (J+i<)  -D’l'l  ( J  +  N)  -  3  .  *3  ( J  +  1  ,  1 )  *C JEE 1  *  A  ( 1 , X+ 1 )  *  Vi  ( 1 ,  .<)  /  A  (  1 ,  1) 

1  +  (3(J  +  1,  X+1  )-3  (o  +  1  ,  l)  *\  (1  ,X  +  1)/A  (1  ,  1)  )  *rt  (  1,u+.<) 

110  CJ  NT IN  JO 

uU  (J)  =302?  1  *C  0  E  F  5  *  D  Y  Y  (J) 

Jii  (J  +  N)  =  J  u  (J+N)  +COE  F 5 
120  CON  i  IN  N a 
F  ZT'J  SM 
EN  3 


S'JJSQJ  TIN  E  PEDZN  V  (I,  l,  ?  NO) 

IMPLICIT  SEAL* 3  (A-H,  Q-2) 

DIMENSION  Y  (3,20)  ,P W ( 4 00)  , A  ( 2 0 ,2 J )  , 3  (20 , 20 ) 


*  *  v  *  ^  *  *  4  *  **  * 


C  5  JbitOU  TINE  CAL-2D  2  Y  3  I7SU3  STORES  IMS  PARTIAL 
C  DERIVATIVES  J?  THE  DIF  F2  RENTAL  L  *J.U  ION  S  PPOVIjSJ  IN 

C  SUBROUTINE  Dlr’JUN.  THE  PARTIAL  DERIVATIVES  (THE  JACoEIAJ) 

C  *Z?E  EVALUAI2J  IN  THE  LAST  SECTION  IN  CHAPTER  5- 


’  *  *  *  *  *  *  j 


\  ************************  ****************  ***++» 


OOMMON/L 1/COE? 1 , CO  EES 
CJ.WON/LJ/N  ,  NT 
CCUON/Li/A  ,  3 
3u  10  A  K  =  1  ,  J 
0  o  4,  b  Ii  -  1  ,  N 

P  i  l  II  +  (  2*  rC  X-  2)  *  .1 )  =COL  7  1*C.jZ?5*S(II  +  1  ,  X  X  +  1  ) 

P..  (11+  ia"‘X.\-1)  *N)  =-CCEF5«j.  *0  (11+  1,  1 )  *COar  1*  \  i  1  ,X.\  +  1,  /A 
?  N  (  II  + 1  *  M  *  S  ♦  ( -  ■*  X  X  -  2  )  *  X  )  =  0  . 

4  J  ?*(  11  +  2  +  14  *  N  ♦  ( 2  *  X  X  -  1  )  =  N  )  =  i  0  (11+  1  *  X  X  +  I  )  -  l,  ( 1 1  +  1 ,  1)  *  A  (  1 ,  .1  ,i  + 
1  *  C  0  a  5 

10  CO  ,i  C  IN  J I 


J  ,  1  i 

•  /  /  •'  ‘  v  '  * 


n  fj 


! 
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susaoun.is  DiCo:i?  (s:j,a,')l, z i) 

DIMENS  ID  M  A  ^20  ,20)  ,UL(2  0, 2  0)  ,  SCALES  (2  0)  ,  IPS  (20) 


:  **  «* «*  *  : 


********  *************.~ 


SU330UTIN  E  CALLED  BY  MATINV  PEfiFCjJS  A  FI3ST  ST  A  0  E  OF 
GAUSSIAN  ELIMINATION  BOUT  IN  E. 


C 


1 

2 

3 

4 

5 


1  0 
1  1 
12 

1  3 

1  4 


1  5 


COMMON  IPS 
J  2=  1 
N=NN 

DO  5  1=1,  N 
IPS  (I)  =  I 

BO  NN  aM  =  0.0 

DO  2  J  =  1  ,N 

UL  (I,  J)  =  A  (I,  J) 

I?  (SO*  NBM-A  ES  (U  L  (I,  J)  )  )  1,2,2 

cOWN&M  =  A3  S  (UL  ( I  , J }  ) 

CO  NT  I N  U E 

IF  (Bu*  OEM)  3,4,3 
SCALES  (I)  =  1.0/BOW  NBM 
GO  TO  j 

CALL  SING  O) 

J2=-  1 

SCALES  (I)  =  0.  G 
CO NT  IN JE 
NHl  a  N- 1 
DO  17  K  =  1  ,  N  Ml 
BIG  =  0.0 
DO  11  I  =  K,N 
IP  =  IPS  (I) 

SIuE  =  A3S  (UL(I?,:<)  )  ^SCALES  (I?) 
IF  (SIDE-BIG)  1  1,  1  1,  10 

Q.X  =  •)!  .ij 

IJXPIV  =  I 

w  0  N  l  I  l1  U  u 

IF  (BIG)  13,12,  13 
CALL  SING  (2) 

J2  =  -  1 

GO  TO  17 

IF  (IDXP  I V  - K )  1  4,  15,  14 

J  =  IPS  (X) 

IPS  (X)  =  IPS  (I D X ? IV  ) 

i.?j  (IJXPIV)  =  J 
K  P  =  IPS  (K) 

PIVJT  =  0 L  ( X ? ,  i\  ) 

X  Pi  =  K*1 
DO  1o  I  =  X  ? 1 , N 
IP  =  IPS  (I) 

EM  =  -JL  {IP,  K) /PIVOT 
JL  (!./,.<)  »  -S3 
DD  1  o  J  a  K ?  1 ,  :1 
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UL(I?,J)  =  *  i.1*'JL  (a?,J) 

1b  CONTINUE 
17  CO  NT  iN  J  z. 

X  ?  =  I  P  S  i  N ) 

:?  (JL(.<?,N))  19,18,19 

16  CALL  S1J 10  (2) 

J2=-  1 
19  PET’J  8N 
END 


SUoKOUTIME  SOLVE  (US  , PUL, 3,  a) 

IMPLICIT  o EAL*8  (A-H,W~E) 

DI.iENSi.3M  PUL  (20,20)  ,  3{2J),a(2J),  IPS  (20) 


'*«**#*****■»*. 


C  SU3SOUTINE  CALLED  31  DIfSUB  PERFORMS  THE  SECOND  STAGS  C 

C  (THE  6  ACri  SUB  S  TIT  JTI  Ctl  PROCESS)  0?  T.iS  GAUSSIAN  C 

C  ELI  HIM  AI  ION  SOUTINE. 
r 


1  **  *c 


C  O.-iJtON  IPS 
N  =  NN 
N?1  =  N+1 
IP  =  IPS  (  1) 

x(i)  *  a  up) 

DO  2  I  «  2 ,  N 
I  P  =  IPS  U) 

I  ;i  1  =  i-i 
SUM  =  0.0 
DO  1  J  =  1,131 

1  SUM  =  SJ.l  ♦  PUL  (IP,S)»X(J) 

2  X(i)  =  o(iP)-su;i 
I?  =  IPS  (N ) 

X(N)  =  X(.«)  /PUL  (IP, N) 

DO  4  I  SACK  =  2,.N 
I  =  NP 1-IBACX 
IP  »  IPS  (I) 

I?1  =  1+1 
S JH  a  0. 0 
DO  J  J  =  I  ?  1  ,  N 

3  SJM  =  GU.1  +  PUL  (IP,  J)  *.<  (U) 

4  X  (I)  =  (X  (*) -SUN) /PUL  (IP,  I) 

?  £  i’U  ON 


SUoSJUIlNS  SINS  (I*HY) 


*«***«  i 


*«*««***i 


i  •  *  a  a  J 


SUBROUTINE  CA-LZ3  JY  D  EC  C  HP  INDICATES  7  HE  EaSOsS  IN 


1 


u  u 


15? 

PSR50J.1ING  TJa  GAUSSIAN  ELI/.  I  NATION  aOUTINI. 


C 


i  «**««  «***4l«« 


1'  FOR  *AI  ( 54  .i  0  H  ATR  IX  WITH  2  c,  HO  ROW  1 N  DECOMPOSE. 

12  r  0  a.1  AT  (  D4  a  3  S  I  NS  UL  A  H  HATRIa  Id  DECOMPOSE.  ZEaO  DIVIDE  Id  _>CL7Z. 

;  Id  IKP5U7.  /.AT5IX  12  dSA5i.iT  G  .t  Li  G  U  L  A  R , 


13 

FOLil  AI 

l54  H 0 

NO  CONVERSE 

NOUT  = 

3 

so  ro 

(1,2, 

J)  , 

iwa  y 

1 

WRITE 

( w  j  j  r 

,11) 

SO  TO 

10 

2 

WRITE 

( d  o  0  T 

,12) 

vjO  r  0 

10 

3 

WRIT  E 

I 

,13) 

10 

RETURN 

u  It  0 

GU3ROU 

TINE 

HAT  I 

dV  (PW  ,  N 

DiHENS 

ION  ? 

a  (400)  ,  A (20 

♦ ** ****** 

***** 

*********** 

3UDRJU7I 

N  £  CALLED 

JY  DIF 

‘  *  *  *  C 


u 

c 


GAUoSIAd  2  L I  H  I  a  AT  I  0  N  Li  GUT  Id  E  OF  THE  HATE  IX  STORED  Id 
THE  ARP. AT  ?W. 

***«**«*.,,> i 


do  i  i  =  i ,  :i 
do  1  J  =  1 ,  :i 

1  a  (i,  j>  =  ii*  {j-  i  )*:;3) 

CALL  DEC*./. P  (drAfULfJI) 

R  ZTU RN 
END 


jUwRJU1I.«L  CD  5  ENT  (  C  0 ,  A  ,  Y ) 
I/.PLICiT  RaAL*3  (A-H,  0-  Z) 

DI/.E  da  ION  A  (20,20)  ,Y  (8  ,  20  J  ,  XS  (2) 


C  SUBROUTINE  CALCULATES  THE  CURRENT  UZNSITi  Bit  El'JAIIud 
C  ( 5o)  . 

v.  *"*  PA  >1  c.  i  .  j  j 

C  w 

c  C3=caa «z;«;  a^p/c  ;<.**<:  : 

C  A-  OliC  jZ i  Li, kZ  lu  COEJFICIINT  *.  AT  3 1 \ ,  A*  C 

u  .  s  C  <j  ^  w  N  **  it  tv  T  ^  0  a  S  A  I*  L  CL  mJl  A  To.  vj  ?  0  L  j  T  3  •  v. 


*********** ^ 


n  r.  n  n  o  n  n 
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c  o  ,1  ,i  .j  /  .t  f  r 

CO.l.i  ON/..b/  XS 

co:ihcn/Lo/c  o^?-. 

C0=-C Jt;4*  (  A  1 1  ,  1)  *X5  (1  )  *A  (  1  , I)  ) 
DU  bo  J  =  1 ,  :i 

CJ=  CD-COE?  ->  *  A  (  1 ,  J  «■ 1 )  *  X  { 1 ,  J) 

50  CONTINUE 
N  ETU  PN 
END 


S  0  3  S  u  U  T I  0  Z  CU'.iCEN  ( C 1  ,  C  2  ,  N  X  >  N  D ,  D  27  A  ,  Z  0  0 1 , 0 1 F 1 ,  i ) 

IMPLICIT  AZAL*3  (A-r.,  O-I) 

01,12 NS -Ob  01  l  60)  ,C2  (60)  ,P.CJT  (2«J)  ,  Jl?1  (20)  ,  Y  (  i,  iO)  ,  15  ,  XI  NO  P  (2U) 

J*#******«*W«***»»»*«*****#***********«********<I**.***  **»*¥***■*¥*  ^ 


C 

C 


5  USr.G  J  Ti  N  Z  Z7ALJA7ZS  TOE  CO  NCZNXuA  TI  ON  AT  ZAC.’i  DZSIAZD 
POINT  IN  X  DISZCTION  E  L  AN  IN  I Z  3  ?0  L  A  T  10  N  SOUTINE  .JOINS 
THE  KN  Oil  ,i  0  J,<  GENERATIONS  AT  TiiZ  COLLOCATION  POINTS. 

- j?  A  r.  A  .1  m  .  R  S  ~“ 

C  1  ,  C2=  CD  NT  E.i  I  LA  TIONS  u?  NITRAIZ  AND  nYDRJXY  lUNJ  , 
RESPECTIVELY*  IN  D I HE  NS  ION L Zoo  UNIT. 

D  E  T  A  =  X  INTZa  i  A  -  ,  DIMENSIONLESS*  0  <C  0  Z  T  A  <  1. 


L 

c 


I  *■*  *  »  : 


CO  H.1  ON /L2/C  OZ  Z  2 
CG.UION/Li/N  *  S  I 
COMauJi/Lb/Xo 
DO  3  0  J=1 , N  X 
S  TA  =D  ZTA  *  (J-  1) 

C A m L  I N  TIT  ( .<  0  *  N T  *  ET  A ,  K  COT  *  D I  r  1  *  \  *.  .i  Z -  ) 
Cl  ( J)  -XiN  TP  (  1)  *XS  (1  )  ♦  X I N  T  ?  (NT)  -1. 

C2  ( J)  =  XiN  TP  (  1)  *  X  5  (2)  +XIN  IP  (NT)  *COZP2 
oo  7o  i , .« 

-  1  (J)  =0  1  (vl)  +  XINT?  (X*  1)  *Y  (1  ,  K) 

7  0  0  ^  (J)  =C2  io)  *-XlNT?  ( X  +  1 )  *i  (1,  h+.i) 

30  Tu  N  I  IN  U  Z 
r  Z  HI  aN 
ZND 


.  li  tt  •*  1 J  .  [  |W  WVU,'  3  f.JCl,  4  f  .U^llD  |  ,*  J.  ,  W,i\  i  | 

1  oT  A  j  .1  A ,  J  o  i  A  o  l  ) 

.OPTIC  IT  j  ( A  -  ii,  0-  E ) 


J'J.  ^  J  -»  O  .  .1«i.  J?  y.  J  .  ^  ^  .  .1  T  |  ,  A  .1  J  .  .1  .  b  ,  .1.  J 

A  .i .  o  .  S  a  ii  .  v.  ,1  A  7  L  HA  .j  jy  0  J  ii.  .0  T  ri  w  i ,,  C:>  .  AJ.  o  -  .  ,i  j 

J?L.  Jy  POT.,  i  • 
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'*******«*  ******* ******* 


u 


3avz5=C3 

IF  (0  3.  i  Z  .  0  »  u  0  1 .  VI  D  .  Z  S  .  l.  T .  J  .  0  1  J ) 
IF  (-  J.  3  E .  J.  U  1  Om  A  N  D  .  Z  5.  LT-  0.  10  0) 
IF  (u5*vjZ«  J-  1  0  0 )  D  u  3  =  0 •  05 
ZS  =  *.S*DZ-, 

FACTO»I=_3/SAV03 

CO  ZF  J  =  0D  if  j/f  ACTOR 

CuEF 4  =  CJ EF4 /FACTO  a 

CO  if  3  =  00  Zf  5 /F  ACTOR/?  AC  TO?. 

dtaj=d:a  j/ i  o j . 

ITAJHI*.  0  0  I  *  J  T  A  0 
D I  AO  rt  A  =  2 . 0  *  0  7  A  J 
J  3'TA rT  =  0 

?s:uaa 

FAD 


Ot, o-  J  .  0  005 
D3S=0 .005 


2  032  30  HO  Z  E.7PVI3  (I,  ROOT,  Jif'l  ,  F  ACT  03,  00) 

ILLICIT  a  Z  Au*  3  (  A  -  H,  0- Z ) 

3 1 0  2  .<  5  i  3  .<  i  (o  ,ZZ)  ,  ROOT  (20)  ,  DIP  1  (00  j  ,  71  NT?  ,20;  ,  01  (zC)  ,  73  (2) 

J*****#««**«»****j)i***J#»»»*^*«**j*«?i******»*!;v**<<!i«*»**»**¥f**^ 


c 


n 


3 0 3300 TI  .<  Z  SV  A- OATES  T HZ  CONCENTRATIONS  AT  TSiZ 
POINTS  IN  lUz,  j3*  COOEDIN  AT  Z  DOZ  TO  THE  CH.VIGZ 
SPLINE  POINT. 


COLiCCAtlON 
OF  THE 


C 

C 


'******«****tt»»*»*«*A******ii  **•*«. 


CO:i:iON/L2/COiF2 
C 0 iU O N / L J / N  ,  NT 

co.ihon/L3/:c  j 

00  7  0  :  =  1  ,  N 
31(1)  =  7  '  1  ,  1) 

7  0  3  7(1  ♦  :i )  =  £  (  1  ,  :  f  N ) 

DO  100  i  =  1 ,  N 
RZ=aOOr  (1+  1)  *  r  ACT  C3 
IF(r.E.GT.I.J)  00  TO  5  0 
C.iL  i  IN  :?t'  (.JD,  NT,  3E,  3 OOT,  DIF  1  ,  TIN  T?) 

7  (1  , 1,  =  ill.  i./  (1)  *XS(1)  +  1 1 ..  T  ?  ( .i  T )  *1.0 
7  ( 1  =  7  I  NT?(  1 }  *X3  (2)  ♦  71  ..TP  (NT)  *CJZF2 

Du  ",  J  r.  =  1  ,  .i 

i  (  1  ,  .  j  =  ;  v  1  , 1 )  ♦TINT?  (  ?.*•  1  )  *  u  7  ( 3 ) 
o  0  7  t  1 , 1  *■  N  J  =  7  { 1 , 1  *■  N)  *  7  I  N  T  ?  \  3  *  1 )  *  3  7  ( X  N  ) 

Ou  . i  100 
9  0  7  ( 1  , 1 )  =  1  . 0 

7  t  1  , 1  N )  =  C  0  i  F  2 
100  CCNCINJi 
I-CIUKN 

i  .*  U 


n  r. 
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F  J  NOTION  A U X l  ( XX , K , J AS  ) 

ILLICIT  5  ZAa*6  (A-H,  t-Z) 

01  IE  NS  ID  N  ti  (/0,2  0 )  ,5  (20,23;  ,  X  X  ( 1 )  ,  U  AS  (20) 


C 

c 


«********«******< 


L. 


C  AUX1  IS  THE  NAJZ  OF  THE  FUNCTION  CALLED  31  ZSYST.1  IN  C 
C  THE  SUufiJUIINa  DIFFUN  TO  FUBwISa  THE  VALUES  OF  EQUATIONS  C 
C  (52)  AND  (53)  .  C 


»♦***«! 


CC.1.1  JN/Ll/CGa?  1  ,CO  EF5 
COMMON/Lz/CO£F2 
CG1.10N/L3/N  ,  NT 
C 0 ilil 0 N / L 4 / A,  3 

C  0  .1  .1  O  N  /  L  3  /  —  u  E  F  3  ,  C  0  N 1  ,  C  0  N  2  ,  3  A  N  A  ,  3  2 1  A 

2AU  (2U)  =-3.  *CQEF1  *  (A  (1  ,  1)  *XX  (1)  *A  (  1,  NT)  *  1 . 0) /A ( 1 , 1 ) -^ ( 1 , Nl ) / A  (  1  ,  1 ; 
1  «C  J  EF  2 

DO  5  00  J  =  1  ,  N 

2Aa  (20)  =  j  AN  (20)  -  3.  ♦COE?  1*  A  ( 1  ,  U  +  1 )  *2  AS  (J)  /  A  ( 1  ,  1)  -  A  ( 1 ,  U  +  1) 

1  **.13  (J«-N)/A{1,1) 

5  00  CONTINUE 

IF  U Aft  (20)  .  LT.  0.  )  £A3(20)=0. 

AUX1  =A  ( 1,  1)  *XX  (1)  *  A  (  1,  NT)  *1  .  ►  (CON2*  (yAE  (20)  **  3 AHA)  -CON  1  * 

1  (IX  (1)  *  *  3  ET  A )  J/COEF3 

DO  t>  00  J=1,j 

600  AUX1=AUX1  *h  (  1,J*  1)  *QA5(J) 

SSTU3N 

END 


FUNCTION  A  J X2 ( XX, K  ,3  AS) 

I  u P L I w  I T  ft  a  A  l.  *  3  ( A—  ri,  a ) 

DINE  NS  UN  A (a0,20)  ,5(20, 20)  ,1(3, 20)  ,XX(1)  ,  LAS  (  1) 


•**  *  *  * «*«* . 


L*t**¥**9*«y*v*«*«*«««**S  »***#***»***.**,#*' 


C 

C  AUX2  IS  THE  U  A  .1  a  OF  THE  FUNCTION  CALLED  o'i  Z  S  Y  Z  I J  I  .< 

C  THa  CAIN  t?SuG3u,l  TO  FU3NI5H  T  rta  VALUES  OF  a<*UATIjNS 

C  (52)  AND  (53)  - 


C 

C 


l***********«*«v***l 


CL  .IN  ON  /L  1/COS?  1  ,CO  EF  5 
C  0 SN 0  N/ L2/C Ua  F2 
C  0,1.1  JN/L  J/N  ,  NT 
C  0  CU  O  N  /  L  4  /  A  ,  u 

Cl,. 1.1  ON  /L  a/  -  0  EF  3  ,  CO  N  1  ,CO  N2  ,  J  ALA  ,  SET  A 
Cj.,.1  DN/LU/i 

Ena(  1)  «  - J.  *C-EF  1*  (A  ( 1 ,  1 )  •!*  (1,  ♦  A  (  1  ,  N  T )  *1.0) /A  (  1  ,  1)  -  A  (  1  ,  NT)  /Ail,1) 
1  *Co  E  F2 

DO  5  00  J  =  1  ,  N 
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3A3  (  1)  =  E  A .»  i  1)  -  3.  *C0E?1*A  (  1 
1  *  i  (1,  J*N)/A  (1,1) 

500  cociruus 

JL  F  Ctl  A5  1 1  ).i.T.  0.)  3.*E(1)=0. 

AUA2=A  i  1 ,  1)  *  a  a  ( 1 )  ♦  A  (  1 ,  ,'J  I)  «1  . 
1  U<  l'»  *  *  3  ET  A )  )/C0  ZF3 

DO  6  00  J  =  1  ,  ,i 

600  AUi2  =  A  Ji2  *A  (  1 ,  J  «•  1)  *’/  (1,0) 
5iT’J3S 
E:i  o 


J*1J  **i1  »  J)/a  (1,  1;  -  A 

( C C :i 2  *  (3  A3  ( 1)  **  Jn.lA)  -CO A  1  * 


IOC 


(j***»*******»aw************4<»***»**«****»******4************«**ai.j 

C  C 

C  THIS  PiiJo^A/  oU  L7  ^3  a'2  U  A  T  10  N  G  {  4  d )  ,  (  4  ■)  )  ,  AND  (52)  i  a  R  0  J*j  H  a 

C  (55)  «xTd  1'tiJ  INITIAL  CONDITIONS  JOINS  SPLINE  COLLOCATION  C 

C  METHOD  FOR  IHE  CASE  CF  THE  ELECTRO!  ITT  CONTA  ININS  CADMIUM  C 

C  IONS.  IN  THIS  CASE,  THE  TWO  SUBROUTINES,  DIF FUN  AND  PEDERV,  C 

C  CALLED  BY  THE  INTEGRATION  SUBROUTINE,  DIFSU8,  HAVE  TC  BE  C 

C  CHANGED.  MEANWHILE,  THE  MAI  N  P20G  RAM  AS  WELL  AS  THE  CTnaR  C 

C  SUBROUTINES  STAY  THE  SAME.  THE  CONCENTRATION  PROFILE  OF  C 

C  THE  CADMIUM  ION  :1AY  BE  OBTAINED  BY  THE  ELECTRONS JT R ALII Y  C 

C  CONDITIO...  THE  CURRENT  DENSITY  DATA  ARE  CALCULATED  IN  IHe,  C 

C  SAME  WAY  EY  EQUATION  (56).  THE  H  E  TE2  UGS  N  EOUS  REACTION  RATE  C 

C  CONSTANTS  OBTAINED  BY  LAST  P3CG2AH  ARE  USED  TO  IDENTIFY 


C  THE  HD  NO SEN  £0  US  REACTION  RATE  CONSTANT  BY  FITTING  THE  C 
C  CALCULATED  CJHcENT  DENSITY  DATA  WITH  THE  EXPERIMENTAL  C 
C  DATA-  C 
C  0 
C  -  ADDITIONAL  NOMENCLATURE  -  C 


C  SM  ALLK  =  HDMDGE  NEOUS  REACTION  RATS  CONSTANT,  1/SEC.  C 

C  D5IGKS  DIMENSION  LESS  HOMOGENEOUS  REACTION  RATE  CONSTANT,  C 

C  05  i  G  K  =  3  Urt  LLX  *  D  MA  X  *  *  2/5.  2b  E  ~5  .  D 

C  D  KS  ?=  D  lil  LNS  10  N  L  ESS  SOLUBILITY  PRODUCT  0?  CADMIUM  HYDROXIDE,  C 
C  DKSP=2.0a-23/C1 B**3.  C 


C**  : 


IMPLICIT  REAL  *8  (A-H,  Q-Z) 

DIMENSION  A  U  0,2  0)  ,S  (2  0,2  0)  ,C  1  (60)  ,C2  ( faO)  ,  C  (2D)  ,  Y  (o  ,20)  , 

1  0  IF  1  (20)  ,DIF2  (2  0)  ,DIF3  (2D)  , ROOT  (20)  , 7  ECT  (20)  , ZnR  OR 

2  YMAX  (20)  ,SAVS  (24, 20)  ,PW  (BOO)  ,  X3  (2)  ,  XX  { 1)  ,  W«t  (3)  ,  BAR 
EXTE  RN  AL  A  J  XL 

COMM  ON/ LI /-UEFI , CO  EF5 
CC  '.MON  /Lk/EOf.F2 
C  •  MO  . .  .  j  ,  N  . 

C.mMu.i/24/h  ,  L 
CCMMUN/L5/XG 


COMMON /Lo/C  QZ  ?  4 
CO. MMON/L  7 /N  GIG,  NEE  ,  E? 
COM.MON/LB/COEF3  ,CO  N1  ,CON2  ,  G  AM  A  ,  BETA 
CGMMON/LR / » 

CUM.ION/L  lO/OBIGK  ,  D  KS  ? 


INPUT  DATA  — 


..  xA  0  (5,310/ 
R  EAD  (5,320) 
R  -  AD  i  5 , 3  3  0) 
R  EAD  (5,340) 


..  D  f  1  ,  0  ,  N  1  .  HU,  3.,  .7  /  .  .4 

v.1  3,  C2E,  DM  AX,  DTIME,  N  (,  EG 
■j  A  M  A  ,  j  t?A  ,  C  GM2  ,GO.i  1  ,  M  M  .A  L  L  K 
N  3  a  vj  ,  N  c*.  ,  E  P 


.  .A  L  C  U  L  .A  T  E  GoME  C  0  G  T  A  N  T . 


:  j 


?j  jG  Sal  - 


COZF1  =  (1.  J02D-5)  /5.  26D-5 
COE  F2  =  C2a/C  1  3 


nr,  fin 
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COZr J=  1 1.  *  JxD-5)  *C1S/D.:A  \/ZZ 
C0£F4=2.  *  1 1.  'i')2D-  5)  *C  1  3*  0  b  43  7_/ J.i  AX /X,  S 

0  3  i.  j  F  =  3 ,,  rt  L  L  *  D.N  A  *  .,.1  A.*.  /  5*  2o  2 
DK3?=(Z.D-2  jj  /(Cl  3**  3.  ) 

DZ=D.1A  a/(Na-1) 

DXAJ  =  (5.  2bi-5)  *DTI  ME/D J AX /DO A a 
TAU1  AX  =  (5.  2  oD-  5)  *TIM  AX/ DM  A  2/D  .1 A ,( 

5  OTA  U=DTAL! 

D  ^  7A  =  1  .  /  v  ii  .v  -  1 ) 

NX=N  *  N  0  +  N  1 
N  ZL  =  2*  J 


INITIAL  DI  il  ZN  SI  0N..Z3S  11.12-0.  - 


7  A  U  =  0. 


-  j  c.'  1 11  .  .  A  —  M  A  a  x  .10  3  rv  N  D  .ix  .<  *  .1 U  .1  D  x  .1  ~  NS  x 0  N  x  —5  j 

11.0  a  IJC3S  MEN!  TO  j;  JSZO  IN  TKZ  INTZG? A  IIOJ 
SUBROUTINE  0  IF  S  03 - 

dxaj.u  =  .  0  0  1 *DTA J 
D T  A 0  ,1i\  -  2  .  0  *  ST  A'J 
*3112(0, JbJ) 

-  £  V  rt  L  u  A  X  2  THE  COLLOCATION  POI.'I  IS  (aCOTS)  OF  I  3 

JACOBI  POLYNOMIAL  OF  OD3Z3  0,  AS  WELL  AS  XH2 
Fid  ST  AND  SECOND  02 hx VAT I V SS  Of  Tr.Z  PoLiNoOI 
A I  r  ri  £  FOOTS. - 

v.  .1  Li*  J  C  %j  ib  x  (  JO,  .i ,  ..  0  ,  1,  rt  3  x ,  0 1  r  1,  D  x  r  2 , 01  *  J  ,  G  .) 

*»NIX2  (b  ,  J  5  0)  (FOOX(J)  ,0=1  , 7 ) 


.CJLATZ  TiiZ  D I SCa ETI2A 71 0 J  COEFFICIZ.iT  MAX 


10=1 

SO  x  J  1=1,..  I 

CA.L  0  :  0  2  3  ( !!  D  ,  N  ,  NO  ,  N  1  ,1  ,  10,0I?1  ,  DI F  x  »  D  x F  .1 ,  EC GT  ,  V  .L 
DO  ID  J  =  1  , N  T 
A(I,J)  =  VwCT  ( J) 

.uLTI  JO  z, 

— —  —  —  —  v^.-ixCox.'i.x  3  ij.i  LL  ^1,  r  Lwij.i  .  1  u  . 


DO  * o  x  =  1  ,  .*  l 

Cd  .  .  x  J  l'  .i  ^  ..  ,  ..v/  1  ,  x  ,  .  J  ,  j  I  i  j.;  3  ,  .i 

Diy  JO  0  =  i  ,  J  . 

DC,J,=7ZCT(J) 


1 ..  lj  J.  r.:Z  V  A  L  0  x  j  jf  A  3  0  i.ix  <Tj  X  j  j  ' 
IN.  _  j  F  \ T  1 0  N  5  J  j  NO  ill  N  X  Ox  .‘003  — — — — 


J  J  X  -> 
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NF=1 

U  S  T  A  nT  =  0 
ISTD ?=0 
HAXD  23  =  7 
DO  47  1= 1 , a  2a 
47  Y.'lAX  (I)  =1.  J 

1  =  0 
n 

C  -  INPUT  THE  INITIAL  CCNUI TIONS - 

C 

DO  50 

i  (1  ,J)  =  U 

50  Y  (1  ,J+N)  =2  0  27  2 
XS  (1  )  =  1. 

XS  (2) =  2DEr  2 
WrilTE  (6  ,  J90  ) 

60  TI2£=TAU*DHAa*DHAX/{5.  2  6D-5) 

v. 

C  - CALCULATE  AMD  PHI  NT  1:12  CL  a  La  NT  DENS  IT  i  AT  ISIS  2  1:12 

C  P  3 1  NT  THE  CONCENTRATIONS  0?  NITRATE  AND  HYDaOXi  IONS 

C  AT  Ill E  NT  COLLOCATION  POINTS - 


CaLL  C u  3  2 N  T  (C  D ,  A  ,  Y ) 

"SITE  (6  ,  ioU)  TIHZ,CD,X3(  1)  ,  XS  (2)  ,  X  F  a  A  G ,  I T  S3  N  ,  23 
NnlTE  (6 ,700j  (  Y  ( 1  ,  L )  ,  L  =  1  ,  N  E  a ) 

C  -  EVALUATE  AND  PAINT  THE  C0NC2N T S A TI C NS  0? 

C  NIl'BATE  AND  HYD30XY  IONS  A I  EACH  POINT 

C  IN  X  DIRECTION  AT  THIS  II H  3 - 

C 

CALL  CO.NCZN  (C  1,C2,  NX  , N  D  , D  ET A  ,  ROOT ,  DIE  1  ,  Y) 
vi SITE  (6 , 2 J 0 )  (Cl  (J)  ,J=  1,  NX) 

<«  a  I T  2  ( o  ,  2  5  0 )  (C  2  ( J)  ,  J  =  1 ,  NX) 

C 

C  - CHECK  IF  THE  I  NTS  S3  ATI  ON  TIHE  AT  THIS  HU. 12  NT  REACHES 

C  THE  HAXllIUa  TINE  TO  32  INTEGRATED - 

C 

IF  (TAU.SE.  TAU11AX)  GO  TO  190 
35  1  =  1+  1 

IF  (I.  EQ.  1 )  GO  TO  9  0 
IF  (ISTOP. £*.  1)  GO  TO  90 
C 


n 

c 

c 


TEST  IF  THE  CO  NCS  NT  3A  TI  ON  0?  NIT3ATE  ION  AT  THE  NTH 
-CuLCCATIOM  POINT, 12,  THE  LAST  ONE  3EF03Z  THE  SPLINE 
PC*  a  I  IS  *  IT  HI  N  J.  OOC  1  Dl.i  ENS  I  JN  LESS  C  D  NCEN  13  a  . I  w  i< 
UNiT  OF  THE  30UNDA3Y  CONDITION-  - 


C 


C 

C 


I?  ((  1.  -Y  ( 1  ,  5)  )  .  LT.  0.  OCO  1)  GO  TO  )0 


IF  THE  TEST  FAILS,  THE  SPLINE  POINT  L3  I  NCI  SAG  2 J , I  2 , 
1  JVi„  FURTHER  INTO  Tri  2  SOLUTION.  - 


» 


c 

c 

c 

c 
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o  A ..  L  w  ri  A  N  vj  i  CO  £..-  d,  r  Av.IOB  ,  0  T  A  U  ,  0  .A  J  HI ,  D  iA  J  H  A  , 

J5TA3T) 


£.7  ALU  ATE  T  iiE  CO  NC  u.<  13  AT  ION  J  J  F  N  IT?  A  72  A  .'I  Zi  HYOFOXT 
IONS  AT  THE  COi.LOCa  II  ON  POINTS  AT  THE  NE'N  CCtJdJIHaTE 
Oil  2  TO  THE  INC  3  EASE  0?  THE  SPLINE  POINT. - 


CALL  EXPAN  0 { Y  , HOOT, DI F 1 , FACTOR , NO) 
C 


INI  EG  3  AT  E  EQUATIONS  ('48)  AND  (4j)  USINS  THE  C  J  NC  E  N  13  AT  1 0  N3 
AT  IhS  COLLOCATION  POINTS  AT  THIS  TINE  TO  OLTAxN  I  .:2 
CONCENTRATIONS  AT  THE  SUCCESSIVE  TINE  - 


90  CALL  DIFSUS  ( N  2  It  ,  TA  U#  Y, S  A7E ,  OTA  J ,0  TAUN I  ,  OTA  UNA  ,  EPS,  NT,  l  H  A  X  , 

1  £.3303,  K?LAG,JSTA3T,NAX0 23  ,  ?  « ) 

I?( XFLAG. LT. 0)  GO  TO  190 
DTA  (JMI  =  . 00 1  * DTAO 
OTA  UN  A  =  2 . 0  *C  TA  rJ 
IF  (I.GE.  2)  Go  TO  93 

IF  (Y  (1  ,  1 )  .  G  1..  OR.  Y  (1  ,  2)  .  ST.  1-. 03.  L  (  1 ,  3)  -  ST.  1.  .O.t.  I  (1  ,  4)  .  ST.  1  . 
1  .03.  DAoS  (  Y  (  1 , 5)  -  1.  )  .  GT.  0.  0004  j  GO  TO  90 

GO  TO  95 
93  N 1=I/IC 
N2=ai*IC 

I?  (I .  E Q.  .12.  0  a.  TA  U .  GE« TA  UNA X)  GO  TO  95 
GO  TO  d 5 


C 


-  CALCULATE  THE  SURFACE  CONCENTRATIONS  - 

95  OLDXS2  =  Xo (2 ) 

X  X  ( 1 J  =  Y  (  1 ,  1} 

ITL3 N=500 

CALL  2SYSTU (aUX2, EP, NSIG, NE'E, XX ,IT23N, AA, 3AR, IER) 

XS(1)  =XX  (1) 

XS  (2)  =-J.  *  CO  EF1*  (A  { 1,  1  )*XX(  1)  *  A  (  1,  NT)  *  1 . 0)  /A  i  1 ,  1 )  -  A  (  1  ,  N  T)  /A  ( 1  , 1 ) 
1  *COEF2 

00  150  J  =  1  ,  N 

XS  (2)  =XS  U)  -  3.  *CCEF1  *A(  1 ,  J*  1)  *Y  (  1 ,  J)  /  A  ( 1  ,  1)  -  A  { 1  ,  J*1) 

1  »Il1,J*S)/l(l,1) 

150  CONTINUE 

IF  (XS  (2)  .  LI.OLOX5  2)  ISTO?=1 
GO  TO  oO 

190  A  3 1 T  E  (a  ,24  0  j  0  ( ,  3  OTAU,  NX,  OTAJ,  DTI  HE,  OH  AX,  AL,SE,  ES,  EPS,  .<  F  La  j, 

1  E? , NSIG 

fl3iiu(tj,245)  SETA,  GA  N  A, CO  Nl, CON  2, SHALL  a 


C 

230 


- TON  .UTS  FOP  INPUT  A  NO  UJTPJT  ST  A  TEN  ENTS - 

FOoH  AT  (/,25  X,  •  NO  3-  ’,  1  X  ,  *C  1  1. 4/3  OX,  90  11 . 4/3  OX,  1  1  .  */jOX,  90  11 
/jOa  ,  90  1  1  .  4/3  0X  ,  9D  1  1  .4/3  OX,  901  1. 4/3  0X,  9*  1.  a/aO  a,  9D  1 1 


J 


* 
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2 

2  40 

1 

2 

3 

4 

5 

2  45 

1 

2 

3 

2  50 

1 

2 

3  10*" 
320 
330 

3  40 
3  50 
360 
3  70 
300 
3  00 
7  00 


/JO  2  ,  00  1  1 .  4/3  Oa,  OD  1  1 . 4/3  02,  301  1 .4/302, 9  0  11.  4/ JO  X,  j'o  11. 4) 
FGSHAi  (//,5X,  'X  INTE3VAL=' , 012- 4,322, ' INITIAL  IAO  STEP  SUE=' 
,012.4/52,  1  MO  Or'  POINTS  I A  2  013=  ',13,232, 

•LAST  TAT  STEP  SUS=  '  ,012.4/52, 'PEAL  Tl;iE  INTERVAL*  ' 
,010.4,242,  'DIFJSION  L  SN  OTH=  '  ,010.  4//  52,  '  A  L  =  *,.-7.2,302 

'  S  li  =  '  ,F7.  2,252  ,'  ZS  =  '  ,  F7 . 4//5  X  ,  '  S?S=  ',010.2,27a, 

'  K  FL  AO=  ',13// 5  2, '2?  =  '  ,  J  1 0 . 2 , 27  X , '  :i  3 10=  ',13) 

FOa.1  AT  (//,  52,  'REACTION  03028  OF  SC3-  ION  =',F7.  2,25a, 

'  4  2  ACT  10  N  0  50  25  OF  OH-  UN  =  '  ,  F  7 . 2//S  2  ,  '  1ST  CONSTANT  =' 

,  0  12.4,252,  '2ND  CONSTANT  =',012.4 
//5  2 , ' REACTION  RATS  CONST AN T=  ',012.4) 

FOSa AT (/, 2 52, 'OH-' ,2X, 9011. 4/302,901 1.4/302,001 1.4/302, 301 1.4 
/JO 2, 9  0  1  1.  4/3  02  , SO  1  1.4/3  0  2, 301  1. 4/3  0  2, 3 0  1  1.  4/ 50 2 , 3 0  1 1 . 4 
/3  02,  4  011.  4/3  OX,  90  1  1  .  4/3  02,  9  01  1.4/30  2,001  1.4/30  2,  ib  11.4) 
FOR  HAT  {413,2/5.  1,  2  E10.  2,1  5,  F3. 4) 


FJ  RH  AT  {  4  0  1  1.4,I7,F7.4) 

FOSSA!  k/F7.  2,3012.  4) 

FO 5.1  AT  (215, 012.4) 

FOBS  AT  (152,  1  OF  10. 6) 

FORI  AT  (//,  uU. 4, 12, D  12.  4,  22, 203  0.10, 215, F2  0.  7) 
POSH  AT (/,  22, qF20.  6/2  X, 6:20. o) 

F0  5H  AT  (22,  '  COLLOCATION  POINTS:  ') 

FOR  .1  AT  (//,  j  a,  'TI.-1E  '  ,  122’  ,  '  C  0  rlaSNT'  ,452,  '  ION  CON 
F osa  AT  (302,  1 0  F 1 0 .  5/302,  10F10.  5) 

STOP 


ENTrATIJN  '  ) 


END 


5  U330U  TINS  01FF0N  (1,72,022) 
iaPLICII  REAl* >3  ( A  -  H ,  Q-  E ) 

01  HE  NS  ION  0  2  2  (20)  ,  YY  (S ,  20)  ,  A  (20, 20)  ,3(20,20),C(20),  22(1)  ,  «A  (3)  , 

1  ^ An  (20) 


C  SUBROUTINE  CALLED  3Y  DIF303  EVALUATES  THE  DERIVATIVES  OF 
C  DEPENDENT  VARIABLES  ST05E0  IN  Y  Y  (1,1)  , 1  =  1,2*  N,  AND  STuSEJ 
C  THE  DERIVATIVES  IN  THE  ARRAY  022.  T  IS  Til 2  INDEPENDENT 
C  VARIABLE. 


EXTERN dL  AO  XI 
COa.10N/L1/CO£F1,COE?5 
Cuaa JN/L2/COEF2 

coaa on /LJ/d , n  t 

c o  aa  ON  /  L 4  /  a  ,  3 
co  an  UN  /*.  7/.«  S.  j  ,  N  r.  C  , :,? 
Coaa JN/LIO/DbIGK,  OKS  ? 
11=2 *N 

JO  3  0  1=  1  , 1 1 
*  An  (I)  =22  (1,1) 

30  CO NT  IN  J  E 
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- CALCULATE  THE  VALUES  C4  THE  SURFACE  [.i=0)  - 

■*X(1)=YY(1,  1) 

I T  E  b  4  =  a  u  0 

u  f\  Ll  uo  /  j  i  it  ( A  U  a  1  ,  -j  ?  /  .i  S  i.  m  i  *  A  ^  ^  A  b  ,  ^  u  b  ) 

CT=-3.  *C0EF  1*  (A  (1  ,  1)  *XX  (1)  *A  (1  ,  .i  r )  *1-  )/A  0 , 1)  -  A  l  1,41)  /A  I  1  ,  1)  *CCEF2 
DO  4 5  1=1,4 

C(I)=u(I*1,1)  *C?  +  B  ( I  ■*■  1  ,N  I)  *C0£  F2 
45  CONTINUE 

- STORE  THE  DERIVATIVES  I A  A  A  4  A  Y  DYY - 

DC  1 20  J=  1  ,  4 

DYY  (J)  =  o  (J»  1,1)  *  X  X  ( 1)  +  d(  J+  1  ,  .,Tj  *  1  . 

DYY  (J*4)  =  C  (J) 

DJ  110  X  =  1  ,  ;; 

D  YY  (J)  =  DYY  (J)  +3  (J  +  1,  <C*1)  *YY  (1,K) 

110  D  YY  (J*N)  =0  YY  (J*  4)  -3  .  *B  { J  *•  1  ,  1)  *CC£F  1  *A  (  1  ,  X  *•  1)  *  Y  Y  (  1 ,  X )  / A  A  1 ,  1) 

1  ♦ ( 3 { J  *  1 , K  *  1 ) -  3 { J ♦ 1 , 1 ) *  A ( 1 , K ♦ 1 ) /A  i  1 ,  1 j ) *  Y  Y ( 1 , ft ♦ A j 

DYY  (J)  =£0£F1  *CCSF5*DYY  (J  ) 

DYY  (J  )  =0  YY  ( J  >  N  )  *C  OE  F  5 

S3?  £3=2.  *  J  DIGS*  (0 .5*  (YY  { 1,  J)  «•  Y  Y  ( 1 ,  J*-  4)  )  -DXSP/Y  Y  ( 1  ,  JO)  /  i  i  .1  ,C  ♦  1.)  ) 
IF  ( S  0  ?  E  A  .  3  1.  0.)  DYY  (J  «■  4)  =  D  Y  Y  (J*-4)  -SUPER 
120  CJSTIUJc 
B  U  .  J  B  4 


S 036 OUT 14  L  PEDES  V  (T, 7,  ?4, 43) 

IMPLICIT  B  E  AL*8  ( A  -  H,  Q-Z  ) 

DIMS 4S IJ4  Y (3, 20)  ,?W  (4  00)  , A  (20,20)  ,3  (2  0,20) 


<■■■*  ********* 


5333CJTI4E  C. 


3  Y  D IF  S  15  STORES  THE  PARTIAL 

PROVIDED 


DERIVATIVES  OF  THE  D  IF  FEW  EMI  AL  EQUATIONS  ?BO  VIDEO  I., 
S13HOUTI4S  j„FF'J4.  THE  PARTIAL  DERIVATIVES  (THE  JACOE 
~  '  . -  THE  LAST  ScCTiON  . . . 


so 

WERE  ZVAu'JAiED 


*.  a  !< ) 


14 


14 


■  **  *  *  ** 


tPTE?  5. 


V. 

C 

c 

'**•#***£ 


COMM  04 / L 1 / O  O  £ F 1 , CO  EF  5 
CUAHOU /LJ/H , H T 
C  O  iL*  C  4  /  .  4  /  (i ,  b 
C  0.1.1  04 /■-  1  0  / 0  a  I  OK  ,  0  XS  ? 

Du  1 0  XK=  1  ,  4 
DO  20  1 1=1,.  i 

J* (11+  (2*.\b-2) *4)  =C0EF 1«CC 2F5*d (II* i ,  XL ♦ 1) 

P  »  l  1 1  ♦  l  2  *  Z  K  -  1 )  *  4  )  =-CC£?5*J  .  *  3  ( 1 1  ♦  1  ,  1 )  *  C  0  E  F  1  *  A  (  1  ,  X  *  I  )  /  A  (  '  ,  1  ; 

IF  l  II.  £*.  is.<)  SO  TC  30 

jU  .0  4  u 

30  £u?=2.*0jISK*(0.5*(Y(1,IIj  *  Y  (  1  ,  1 1  ♦  4  )  )  -  D  ,\  „  ?/  Y  ^  1  ,  1 1 »  .i ;  /  i  k  1  ,  1 «.  *•  4  ;  j 

IF(JU?.UI.  0.  )  P  *  ( II  ♦  (  2  ♦K  K-  1)  *:.)  =2*  ( 1 1  *■  (  Z  *  X  K  -  1 )  *4  )  -  DEIS''. 

PA  (  II  *2  •ft  *44  ( 2  *  K  X  -  2  )  *4)  =0. 


40 


1 
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—  *  A  (  1 ,  X  K  *  1 )  /A  1 1  ,  1 )  ) 

♦C0SF5 

Ir(II.U.dH)  JO  TC  u  0 
JO  TO  20 

60  IF(  jjp.  ji.  o- )  ?w  (n*2*M*.s»  (2*k;\-  ij  *:i)  =?«  (  *:;♦  (2*aK-  i )  «..)  -2.  * 

1  D3ISX*  (.  6  *-2.  (  1,  11+:,)  / 
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